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T  temperature  (K) 
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a'  ratio  of  outer  to  inner  magnetic  coil  radius 

0  ratio  of  plasma  kinetic  pressure  to  magnetic  pres¬ 
sure 

/?'  quotient  of  coil  length  and  twice  the  inner  coil 
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7  power  fraction 
A  difference,  change 
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T)  efficiency 

A  volume  fraction  of  the  magnet  devoted  to  produc¬ 
ing  the  field  (around  90%);  Bremsstrahlung  wave¬ 
length  (A) 

p0  permeability  of  free  space  (po  =  4x  x  10-7  H/m) 
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1  INTRODUCTION 


The  demands  placed  on  space  transportation  systems  for  reduced  mission  time  and  increased  payload  capacity 
require  the  development  of  alternate  and  advanced  space  propulsion  concepts.  Deep  space  missions  will 
require  high  specific  impulse  ( I,p )  and  large  changes  in  velocity  (AV).  The  compact  energy  storage  of 
nuclear  fis^.on  systems  and  advanced  fusion  reactor  concepts  make  them  ideal  power  source  candidates  for 
future  space  applications.  Fusion  is  particularly  attractive  because  it  produces  more  energy  than  fission, 
and  does  not  produce  long-lived  radioactive  waste.  Table  1  shows  some  possible  energy  sources  for  space 
propulsion.  With  the  exception  of  antimatter  annihiliation,  fusion  has  the  highest  specific  energy  of  any 
energy  source. 


Table  1:  Comparison  of  Energy  Sources 


Energy  Source 

Specific  Energy 
[J/kg] 

Chemical  (O2/H2) 

1.5  x  107 

Chemical  (02/Be) 

2.6  x  107 

Atomic  Hydrogen  (H+H=H2) 

2.2  x  108 

Metastable  Helium  (He*) 

4.6  x  108 

Nuclear  Fission  (100%) 

8.0  x  1013 

Nuclear  Fusion  (D-3He) 

3.5  x  1014 

Antimatter 

9.2  x  1016 

The  purpose  of  this  study  was  to  produce  a  FORTRAN  computer  code  to  analyze  the  performance  of 
fusion  reactors  as  propulsion  concepts.  The  code  simplifies  the  task  of  comparing  thruster  performance 
while  varying  individual  or  multiple  reactor  parameters.  Shown  in  Fig.  1,  a  generalized  fusion  reactor 
system  was  modeled,  establishing  the  power  balance,  fuel  flow,  power  input,  and  reactor  breakeven  and 
ignition  conditions. 

Two  potential  candidate  fusion  reactor  concepts  were  modeled  and  inserted  into  the  reactor  system 
model:  (1)  the  Dense  Plasma  Focus  (DPF);  and  (2)  the  Field-Reversed  Configuration  (FRC).  The  user 
selects  one  of  the  reactor  types  and  the  fusion  fuel  combination.  The  modular  design  of  the  program  allows 
the  user  to  add  a  subroutine  with  the  model  of  a  different  reactor,  if  desired.  For  a  fusion  thermal  rocket, 
the  propulsion  relations  were  established  between  fusion  output  and  rocket  performance,  system  mass,  and 
mission  capability.  By  varying  fusion  plasma  temperature,  fuel  mixture  ratio,  and  mission  AV,  performance 
values  for  comparision  are  generated,  including  thrust,  jet  power,  jet  specific  power  (a),  thrust-to- 
weight  ratio,  and  payload-to-system  initial  mass  ratio.  The  output  file  can  be  analyzed  directly  with  text 
explanations  or  numeric  columns. 


1.1  Fusion  Fundamentals 

Fusion  is  the  process  that  powers  the  stars.  Fusion  occurs  when  two  light  atomic  nuclei  join  to  form  a  heavier 
nucleus  and  a  lighter  nucleus  with  high  kinetic  energy.  To  achieve  fusion,  the  distance  between  nuclei  must 
be  small  enough  for  the  short-range,  strong  nuclear  attractive  force  to  exceed  the  long-range  Coulombic 
repulsive  force.  Fig.2  [1]  shows  the  Coulomb  potential  barrier  as  a  function  of  distance  from  the  nuclear 
core,  where  Z\  and  Zi  are  the  atomic  numbers  of  fuel  elements  1  and  2  respectively,  e  is  the  charge  of  an 
electron,  and  Ro  is  the  distance  between  the  centers  of  the  nuclei  where  the  strong  nuclear  force  becomes 
dominant.  Fusion  is  easier  with  light  elements  because  the  nuclei  have  fewer  protons,  and  therefore  the 
Coulomb  repulsion  force  is  less.  The  trick  then  is  to  bring  nuclei  together  with  sufficient  energy  to  overcome 
the  initial  repulsion. 

Fusion  reactions  result  in  heavier  nuclear  products  that  are  more  stable  than  the  light  atomic  nuclei 
reactants,  because  their  binding  energy  is  greater.  Fig.3  [2]  shows  the  change  in  binding  energy  per  unit 
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Figure  1:  Fusion  Rocket  Schematic. 


Figure  2:  Coulomb  Potential  Barrier 


atomic  mass  as  a  function  of  nuclear  atomic  mass.  When  two  light  elements  fuse  to  form  a  heavier  element, 
the  the  binding  energy  per  unit  atomic  mass  increases,  leading  to  an  increase  in  stability. 

Fusion  products  have  a  combined  rest-mass  lower  than  the  sum  of  their  initial  rest-masses.  This  mass 
difference  is  converted  to  kinetic  energy  of  the  products  in  accordance  with  the  equation  Ej  =  A  me2  [3], 
where  c  =  3  x  10s  m/s  is  the  speed  of  light. 

1.1.1  Fusion  Fuels. 

In  a  gas,  fusion  does  not  necessarily  occur  during  each  particle  collision.  The  probability  of  a  fusion  reaction 
between  two  particles  can  be  described  in  terms  of  a  cross  section,  <r,  of  one  fuel  particle  to  the  other.  The 
cross  section  for  particle  systems  in  which  one  fuel  element  is  stationary  while  the  other  impinges  (target 
system)  is  different  than  the  cross  section  for  particle  systems  in  which  both  fuel  elements  have  the  same 
average  kinetic  energy.[4]  The  cross  section  also  depends  on  the  mass  of  the  particles,  their  atomic  charge,  the 
nuclear  neutron-to-proton  ratio,  and  the  nuclear  spin,  but  generally  is  highest  for  particles  with  low  atomic 
charge  and  high  energy.  In  choosing  a  fusion  fuel,  one  must  not  only  consider  the  engineering  difficulty  of 
fusing  that  fuel,  but  factors  such  as  the  fusion  reaction  products,  materials  limitations,  and  fuel  availability. 

The  most  commonly  examined  fusion  reactions  for  light  elements  are  shown  in  Table  2  [4, 5, 6, 7].  The 
fusion  fuel  with  the  highest  cross  section  at  the  lowest  energy  is  an  equal  mix  of  deuterium  and  tritium.  The 
temperatures  required  for  significant  rates  of  D-T  fusion  peak  around  5  keV  [6],  or  roughly  fifty  eight  million 
degrees  Kelvin.  D-T  fusion  is  the  easiest  to  achieve,  but  the  reaction  gives  80%  of  its  energy  to  neutrons, 
which  carry  no  charge. [5]  While  adequate  for  ground-based  power  generation,  the  massive  shielding  necessary 
to  protect  the  system  components  makes  D-T  a  poor  fuel  choice  for  space  applications. 

Pure  deuterium  fusion  (D-D)  requires  higher  temperatures  to  achieve  significant  fusion  rates  and  is 
therefore  technologically  more  difficult,  but  as  a  benefit  it  releases  only  34%  of  its  energy  to  neutrons.  In 
addition,  a  virtually  inexhaustible  supply  of  deuterium  can  be  extracted  from  normal  seawater,  and  the  D-D 
reaction  does  not  require  the  handling  of  radioactive  tritium.  Shown  in  Table  2  sue  the  two  branches  of  D-D 
fusion,  which  occur  with  almost  equal  probability.  Although  tritium  is  one  of  the  products  of  D-D  fusion, 
it  can  either  be  treated  as  waste  or  burned  with  deuterium  (“catalyzed”  D-D).  More  advanced  reactions, 
such  as  p-8Li  and  p-uB  produce  no  neutrons,  but  require  temperatures  as  high  as  kT  =  100  keV  to  produce 
significant  fusion  rates. [8] 

D-3He  appears  to  be  the  most  promising  fusion  fuel  choice  for  propulsion.  [9,10,11,12,13,14,15,16]  The 
reaction  products  are  protons  and  alphas,  and  at  kT  as  30  keV,  the  heating  from  fusion  products  equals  the 
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Table  2:  Standard  (D-T)  and  Advanced  Fusion  Reactions 


Name 

Fusion  Reaction  (MeV); 
Abbreviated  Form 

Energy  Release  E/jk 
(MeV;J) 

D-T 

D+T  —  o(3.5)  +  n(14.1) 
T(d,n)4He 

17.6;  2.82  x  10~13 

D-D(n) 

D+D  -3He(0.82)  +  a(2.45) 
D(d,n)3He 

3.27;  5.24  x  10~13 

D-D(p) 

D+D  —  T(1.01)  +  p(3.02) 
D(d,p)T 

4.03;  6.46  x  10~13 

D-3He 

D+3He  — *■  a(3.87)  +  p(  14.68) 
3He(d,p)4He 

18.35;  2.94  x  10" 12 

p-8Li 

p+8Li  -3He(2.3)  +  a(1.7) 
8Li(p,a)3He 

4.0;  6.41  x  10" 13 

p-nB 

p+nB  — »  3  a(2.88  ea.) 
llB(p,3a)4He 

8.64;  1.38  x  10~1J 
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radiative  heat  losses.  Protons  and  alphas  are  charged  particles,  so  their  motion  can  be  directed  magnetically. 
With  deuterium  and  helium-3  fusion,  the  only  energy  released  to  neutrons  occurs  from  side,  “parasitic” 
D-D  reactions,  and  it  may  be  possible  to  suppress  the  parasitic  D-D  and/or  enhance  the  D-3He  fusion  cross 
sections  by  aligning  the  spins  of  the  nuclei[17,18,19,20,21].  Once  D-3He  fusion  has  been  achieved,  the  only 
remaining  issue  will  be  availability  of  3He.  Currently  this  isotope  of  helium  is  rare  on  earth,  and  the  only 
source  occurs  through  the  beta-decay  of  tritium.  However,  it  appears  that  an  economically  feasible  source 
of  3He  exists  in  the  soil  on  the  Moon’s  surface.  [22,23,24,25]  Mining  of  lunar  3He  might  help  spur  along  the 
development  of  space  assets  and  controlled  fusion  power. 

1.1.2  Reaction  Rates  and  Power. 

A  plasma  is  a  mixture  of  charged  and  neutral  particles  exhibiting  collective  behavior.[26]  Plasmas  are  created 
by  heating  gasses  to  high  temperatures  such  that  gas  molecules  begin  to  dissociate  and  ionize.  In  a  plasma 
with  fuels  j  and  k  at  ion  densities  n;  and  n*,  and  at  energy  kT,  the  “reaction  rate”  depends  on  the  density  of 
each  element  of  the  fuel  and  the  cross  section  for  the  j,  k  combination.  The  average  product  of  cross  section 
and  particle  velocity,  <<rv>,  is  specific  for  a  given  fuel  at  a  given  temperature.  The  fusion  reaction  rate  is 
determined  by  [8] 


RRj.k  — 

1  +  6j,k 

(i) 

6j.k  = 

r  o  j*k 

1 1  ;=* 

(2) 

Tables  of  <av>  values  for  particular  reactions  over  a  plasma  temperature  range  are  available. [4, 5, 6], 
and  Fig.  4  shows  <<rv>  for  some  of  the  reactions  listed  in  Table  2.  When  the  energy  Ejjh,  released  by  a 
particular  fusion  fuel  reaction,  is  multiplied  by  the  reaction  rate,  one  has  the  power  density  of  the  reaction 
which  is  a  convenient  figure  for  evaluating  the  usefulness  of  fusion  fuels.  When  multiplied  by  the  plasma 
volume,  Vp,  the  power  density  gives  total  power  produced  by  the  fusion  process.  The  fusion  power  is: 

Pf  =  RRj'kVpE/jii  (3) 


1.1.3  Radiation. 

The  two  dominant  types  of  radiation  produced  in  a  magnetically  confined  fusion  plasma  are  bremsstrahlung 
and  synchrotron.  [3]  Bremsstrahlung,  which  means  “braking  radiation”  in  German,  occurs  when  charged 
particles  collide  with  other  charged  particles  and  accelerate,  emitting  photons  in  the  process.  Synchrotron 
radiation  is  due  to  the  gyration  of  electrons  around  magnetic  field  lines  The  total  radiation  power  can  be 
expressed  [8]  in  terms  of  plasma  temperature  and  magnetic  pressure  parameters  as 

Pir  +  Pcy  =  (cy/kT,  +  drL——(kTt)2K^j  n}Vp  [W]  (4) 

where 

c  =  5.35  x  KT31/*  kiZi 

„  =  MxlO-"(l  +  Jgr)(l  +  f£) 

f*  = 

j 

Hi 

n. 


(5) 

(6) 

(7) 
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i%i  [ions/cm-3]  is  the  ion  density  in  the  plasma;  kTt  and  kTi  [keV]  are  the  respective  electron  and  ion 
energies  in  the  plasma;  /?  is  the  ratio  of  plasma  pressure  to  magnetic  pressure;  f,  represents  the  number  of 
electrons  present  due  to  ionization  of  different  elements,  where  kj  is  the  fractional  density  of  fuel  element  j. 
Ke  describes  the  fraction  of  synchrotron  radiation  absorbed  in  the  chamber  walls,  (it  depends  on  “plasma 
depth,”  involving  system  geometry,  plasma  electrical  resistivity,  reflectivity  of  the  chamber  walls,  and  plasma 
/?.[8])  In  high  temperature  fusion  systems,  synchrotron  emmision  is  very  large,  and  operation  is  not  possible 
unless  Kc  is  kept  below  about  ten  percent. 


1.1.4  Magnetic  Confinement 


There  are  many  different  confinement  schemes  being  developed  to  harness  the  energy  of  fusion.  The  three 
main  types  being  studied  are:  magnetic  (MCF);  inertial  (ICF);  and  electro-static.  The  Field  Reversed 
Configuration  and  Dense  Plasma  Focus  reactor  designs  considered  as  test  cases  are  magnetic  confinement 
schemes. 

Controlled  fusion  is  difficult  to  achieve  because  it  requires  maintaining  very  high  temperatures  and 
pressures.  To  achieve  fusion  by  heating  and  containing  the  fuel,  the  plasma  particles  must  have  limited 
contact  with  the  reactor’s  walls.  Otherwise,  the  plasma  is  cooled  when  it  transfers  its  kinetic  energy  to  the 
walls. 


[7]: 


The  kinetic  pressure  of  a  plasma  is  the  sum  of  the  kinetic  pressures  of  the  charged  particle  components 


P  -  p>  =  n*kT*  +  nikT'  +  nakT<*  •  • 


(8) 


At  equilibrium,  the  plasma  kinetic  pressure  must  equal  the  confining  magnetic  pressure,  defined  as  B3/2fiQ.[7] 
For  MCF,  the  magnetic  field  determines  the  plasma  pressure,  and  thus  the  plasma  density.  The  fusion 
reaction  rate  depends  on  the  strength  of  the  magnetic  field  indirectly,  through  the  particle  density  and 
plasma  temperature. 


1.1.5  Ionization 


All  the  particles  leaving  the  fusion  plasma  have  energies  in  the  keV  to  MeV  range,  so  it  is  very  likely  that 
they  are  completely  ionized.  The  energies  of  ionization  for  the  propellant  can  be  calculated  using  the  Bohr 
model  for  excitation  [2], 


U .  - 

*  "  (47rco)22AJn2 

n 


(9) 


where  Z  is  the  propellant  atomic  charge,  and  n  is  the  state  of  each  electron  to  be  ionized.  Monatomic 
hydrogen  requires  13.6  eV  to  be  ionized,  so  at  keV  and  MeV  temperatures  hydrogen  ionization  will  be 
virtually  complete.  However,  if  the  hydrogen  temperature  drops,  complete  exhaust  ionization  becomes  less 
certain. 

Using  a  form  of  the  Saha  Equation  for  hydrogen  [7],  we  can  find  the  degree  of  ionization  at  lower 
temperatures.  Setting  the  total  heavy  particle  density  equal  to  the  ionized  density  plus  the  neutral  particle 
density  (n«  =  +  n„),  we  can  find  the  degree  of  ionization  (m/nt)  to  be 


Hi  _  (1  +  4z)1/2  -  1 
n,  —  2x 


(10) 


where  x  is  defined  as 

nth3  exp(Ui/kT) 

X  ~  (2*me*r)3/2 


(11) 
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and  h  is  Planck’s  constant,  m,  is  the  electron  mass,  I/,-  is  the  ionization  potential  (13.6  eV  for  hydrogen),  k 
is  Boltzmann’s  constant,  and  T  is  the  temperature.  If  kT  is  expressed  in  eV,  then  we  can  simplify  x  to 

x  =  3.313  x  10"Mnt(/fer)-3/3  exp {Ui/kT)  (12) 

Fig.  5  shows  the  degree  of  ionization  versus  plasma  temperature.  For  nt  =  1015  cm-3,  ionization  is 
essentially  complete  at  a  few  eV. 


Ion  density 
(cm-3) 

o  1  el  3 
•  1  el  5 
▼  1  el  7 
»  1  el  9 
a  1  e2 1 
■  1  e23 
a  1  e25 


Figure  5:  Ionized  Fraction  for  Hydrogen,  Based  on  Saha  Equation. 


1.2  Rocket  Propulsion  Fundamentals 

A  rocket  provides  thrust  by  Newton’s  Second  and  Third  Laws  and  the  principle  of  conservation  of  momentum. 
Newton’s  Second  Law  states  that  the  acceleration,  a,  of  a  body  with  mass,  m,  is  due  to  an  applied  force 
F  (F  =  ma).  Newton’s  Third  Law  states  that  every  action  must  have  an  equal  and  opposite  reaction.  If 
propellant  is  accelerated  in  one  direction  by  a  propulsive  force,  the  rocket  feels  the  same  force  or  thrust 
in  the  opposite  direction.  The  rocket  ejects  propellant  in  the  opposite  direction  of  desired  motion,  and 
the  momentum  imparted  to  the  propellant  must  be  equal  in  magnitude  but  opposite  in  direction  to  the 
momentum  imparted  to  the  rocket,  such  that  total  momentum  is  conserved. 

1.2.1  Rocket  Equation 

The  change  in  velocity  of  a  body  is  the  integral  of  the  body’s  acceleration  over  time.  If  a  rocket  of  mass 
m  exhausts  propellant  of  mass  dm  at  velocity  ve,  the  rocket  velocity  increases  by  dv.  By  conservation  of 
momentum, 

mdv  =  —vedm  (13) 

Integrating  Eq.  (13)  and  raising  both  sides  to  powers  of  e  gives  the  well  known  Rocket  Equation  [27] 


where  the  mo  and  mj  are  the  initial  and  final  rocket  masses,  respectively. 

The  specific  impulse  (Itp  =  vt/g)  of  the  rocket  measures  how  efficiently  the  rocket  uses  its  propellant  to 
produce  thrust.  One  wants  the  I,p  to  be  as  high  as  possible  since  any  extra  mass  carried  aboard  the  rocket 
decreases  the  amount  of  payload  allowed. 
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The  Rocket  Equation  is  used  to  determine  how  much  of  the  rocket’s  total  mass  can  be  payload.  It  is  most 
desirable  to  have  the  ratio  mo/m/  as  close  as  possible  to  one,  meaning  the  greatest  possible  fraction  of  the 
rocket  is  payload.  Av  is  essentially  Axed  for  a  given  mission.  Therefore,  we  need  to  make  the  denominator 
in  the  exponent  of  Eq.  (14)  to  be  large  compared  to  the  Av,  i.e.  have  a  high  I,p.  A  propellant  with  low 
molecular  weight  is  easier  to  accelerate  to  high  exhaust  velocities,  thus  improving  the  payload  capabilities 
of  the  rocket. [28] 


1.2.2  Nozzle  Theory 

The  objective  of  a  nozzle  is  to  convert  random  thermal  energy  into  directed  kinetic  energy.  Knowing  char¬ 
acteristics  of  the  nozzle,  such  as  exit-to-throat  area  ratio  and  length,  allows  us  to  describe  thrust-generating 
ability. 

There  are  several  assumptions  customarily  made  in  deriving  theoretical  nozzle  performance:  (a)  working 
fluid  is  homogeneous  and  invariant;  (b)  working  fluid  obeys  ideal  gas  laws;  (c)  there  is  no  friction;  (d)  flow  is 
adiabatic;  (e)  propellant  flow  is  steady  and  constant;  (f)  no  reactions  are  occurring  in  propellant  to  change 
its  composition  or  enthalpy;  (g)  all  velocities  exiting  are  axially  directed,  allowing  one-dimensional  analysis; 
(h)  gas/plasma  velocity,  pressure,  or  density  is  uniform  for  any  nozzle  cross  section.[29]  These  assumptions 
sure  idealizations  of  nozzle  conditions,  but  are  usually  accurate  to  ten  percent  or  less. 

Formulas  to  be  used  include  conservation  of  energy  which,  under  the  above  assumptions,  reduces  to: 


[‘41.44] 

The  ideal  gas  law:  p  =  nkT ;  The  isentropic  flow  process: 

Msf  40 


> 


(15) 


(16) 


where  7  is  the  ratio  of  the  specific  heat  for  constant  pressure  Cp  and  the  specific  heat  for  constant  volume 
Cv ;  and  the  acoustic  velocity  definition 

a=v^T  (17) 

Using  the  assumptions  and  formulae  listed  above,  the  velocity  at  any  point  in  the  nozzlp  can  be  found  in 
terms  of  the  enthalpy  difference  of  another  point: 


vj  =  ^2(hi  -  h\)  +  vl 


N 


1-(E)  ’ 

7-1  \PoJ 


(18) 


In  this  case,  point  number  1  has  been  chosen  to  be  the  inlet  stagnation  pressure  and  temperature,  po  and 
To-  As  can  be  seen,  the  exhaust  velocity  i>2  reaches  a  maximum  when  the  pressure  ratio  ps/po  goes  to  zero, 
corresponding  to  an  infinite  expansion  of  the  nozzle. 

If  a  converging/di  verging  nozzle  is  employed,  flow  at  the  “throat”  is  “choked” ,  indicated  with  a  superscript 
*.  The  flow  velocity  at  the  throat,  v*,  is  equal  to  the  local  acoustic  velocity  a.  the  pressure  and  temperature 
at  the  throat  correspond  to  the  maximum  isentropic  mass  flow.  If  the  pressure  ratio  after  the  throat  is 
decreased  by  expanding  the  nozzle,  the  propellant  velocity  exiting  the  nozzle,  vtx  will  exceed  v*,  producing 
supersonic  flow. 

Thus,  including  each  type  of  particle  exiting  the  nozzle,  the  thrust  of  the  nozzle  is 

F  =  ^miVtx<i  (19) 
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1.2.3  Rocket  Performance 


The  total  system  mass  will  consist  of:  the  reactor;  necessary  shielding;  propellant,  fuel,  structures  and  tanks; 
the  radiator;  the  nozzle  (magnetic  or  conventional);  and  the  payload.  Since  everything  but  the  payload 
is  “overhead,”  it  is  desirable  to  reduce  the  system  mass  as  much  as  practical  to  increase  payload  and/or 
decrease  mission  time. 

Given  the  thrust,  specific  impulse,  and  initial  mass,  it  is  possible  to  determine  the  performance  of  the 
rocket.  The  thrust-  to-weight  ratio  is  F/(mog)-  The  jet  specific  power,  a,  is  defined  as 

mp  _  mp 
“  ”  Pjt,  ~  *FI. p 


We  define  the  starting  mass,  mp,  to  be  the  sum  of  masses  of  the  ship,  m,  (which  includes  the  payload, 
reactor,  power  system,  etc.),  the  propellant,  m?rop,  and  the  propellant  tank,  mtank- 


mo 

—  mt  *+■  TYlprop  ■+*  TTltank 

(21) 

m, 

— —  771 0  T7 Iprop 

(22) 

prop 

=  Wlproptfire 

(23) 

For  missions  with  long  firing  times  and  high  propellant  flow  rates,  the  propellant  tank  mass  can  be  large 
compared  to  the  ship  mass.  If  we  define  the  tank  fraction  as  Xtonk  =  mtankfmprop  and  insert  Eqs.  (21), 
(22),  and  (24)  into  the  Rocket  Equation,  Eq.  (14),  we  can  solve  for  the  mission  firing  time  as  a  function  of 
A v: 

m.  (e^/«.  -  l) 


</,r*  thprop  1  -  Xtank(eAvtv‘  -  1) 

Given  the  payload  mass,  some  fraction  of  the  ship  mass,  the  payload  mass  fraction  at  launch  is  simply 
mpayload/mQ- 


(24) 


1.2.4  Magnetic  Nozzle 

Exhaust  gases  from  a  fusion  rocket  are  likely  to  be  partially  or  completely  ionized,  thus  a  magnetic  nozzle[30,31,32] 
is  an  logical  choice.  Fig.  6  shows  the  magnetic  field  required  for  thrust  conversion  of  different  exhaust  den¬ 
sities,  as  a  function  of  temperature. 

An  azimuthal  magnetic  coil  is  one  candidate  for  a  nozzle  substitute.  However,  there  exist  some  engineering 
issues  that  must  be  resolved  before  operational  used  of  magnetic  nozzles  is  assured.  [31] 

1.  At  plasma  temperatures  below  1  eV,  the  plasma  /?  must  be  below  unity  to  prevent  unacceptable 
amounts  of  resistive  cross  field  mass  transport  (particles  not  following  field  lines);  at  10  eV  and  /?=  1, 
this  is  not  a  problem. 

2.  In  the  temperature  and  particle  density  range  we  wil!  most  likely  be  operating  (reactor-like  conditions), 
nozzle  generation  of  radiation  cm  be  severe  [31]  (on  the  order  of  a  few  MW/m2  for  a  1  m  diameter 
nozzle  throat  and  densities  up  to  1018  cm-3)  and  have  an  prohibitive  impact  on  wall-loading  limitations 
and  cooling  requirements. 

3.  Due  to  the  shape  of  the  magnetic  field  lines,  a  “detachment”  problem  exists  whereby  particles  follow 
the  magnetic  field  lines  around  and  deposit  their  energy  back  on  the  vehicle  (negating  any  thrust 
benefit).  One  of  the  ways  to  counter  the  detachment  problem  is  for  the  plasma  entering  the  nozzle  to 
have  a  field-free  core,  or  similarly  a  region  in  which  the  core  field  lines  doubles  back  on  themselves. 

Fig.  7  shows  characteristics  of  a  Meridional  Magnetic  Nozzle. [31] 
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Figure  6:  Magnetic  Nozzle  Field  vs.  Exhaust  Temperature 
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Figure  7:  Meridional  Magnetic  Nozzle  Characteristics 
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2  FUSION  PROPULSION  SYSTEM  MODEL 

2.1  Fusion  Rocket 


In  the  thruster  design  considered  for  this  study,  the  propellant  and  the  fuel  are  separate,  so  the  lightest 
possible  propellant,  hydrogen,  may  be  used.  The  high  temperature  plasma  from  the  reactor  is  used  to 
dissociate  and  ionize  the  diatomic  hydrogen  into  protons  and  electrons.  The  degree  of  propellant  ionization 
will  depend  on  the  final  temperature  and  pressure  of  the  propellant.  If  a  large  fraction  of  the  propellant  is 
ionized  and  the  fusion  fuel  chosen  such  that  the  reaction  products  are  charged  particles,  we  can  employ  a 
magnetic  nozzle  to  direct  the  charged  particles  and  increase  thrust. 

Rocket  systems  employing  separate  propellant  and  fuels  are  usually  limited  by  specific  power  (a)  instead 
of  specific  impulse  (f«P)[33j  Specific  impulse-limited  systems,  such  as  conventional  chemical  rockets,  are 
typically  high-thrust  vehicles,  and  since  the  fuel  reaction  products  are  the  propellant,  the  amount  of  reaction 
energy  imparted  to  each  particle  is  fixed.  With  a  specific  power-limited  system,  more  power  may  be  added 
to  the  propellant  particles  until  the  power  source  becomes  too  heavy.  Fusion  promises  to  be  an  attractive 
propulsion  power  source  because  even  small  reactors  should  produce  large  powers. 

A  mathematical  model  of  a  fusion  propulsion  system  is  necessary  to  see  the  effects  of  changing  one 
parameter  versus  another.  The  system  model  describes  the  use  of  the  reactor’s  power,  how  the  radiation 
and  thermal  heat  is  dissipated,  method  of  electrical  power  generation,  shielding  requirements,  thrust  pro¬ 
duction,  system  mass,  and  rocket  performance.  It  allows  performance  predictions  and  system  performance 
optimizations. 


2.2  System  Description 

Fig.  1  shows  the  component  diagram  of  the  fusion  propulsion  system  considered  in  this  study.  The  heart  of 
the  fusion  propulsion  system  is  the  reactor,  where  the  fusion  fuel  is  burned  and  high  temperature  charged 
particles  are  produced.  The  reactor  also  produces  electromagnetic  radiation  and  neutrons  that  heat  the  walls 
and  other  components.  Cold  propellant  cools  these  components  by  heat-exchange,  then  mixes  with  charged 
particles  escaping  from  the  fusion  plasma.  The  heated  propellant  is  expanded  through  a  nozzle  (conventional 
or  magnetic,  depending  on  the  propellant  temperature),  providing  thrust  for  the  vehicle. 

If  neutronic  reactions  occur  during  the  fusion  process,  shielding  will  be  necessary  to  protect  sensitive 
components  and/or  the  payload.  As  the  number  of  fusion  neutrons  produced  is  reduced,  the  mass  of  the 
shielding  can  be  similarly  reduced. 

2.2.1  System  Power  Balance 

The  fusion  propulsion  system  may  be  modeled  in  terms  of  its  energy  flow  rates.  A  schematic  of  the  system 
power  balance  is  shown  in  Fig.  8.  [8]  Fusion  power  Pj  is  released  by  the  fusion  reaction  within  the  plasma. 
The  fusion  power  is  distributed  between  charged  particles  and  neutrons,  Pe  and  Pn,  respectively.  A  fraction, 
y,  of  charged  particle  power,  P«,  is  deposited  in  the  plasma,  heating  the  particles  and  electrons.  If  injected 
power  Pinj  is  needed  to  sustain  the  reaction,  the  total  output  power  is  the  sum  of  the  fusion  power,  Pj ,  plus 
the  injected  power,  Pinj- 

The  fusion  plasma  is  cooled  by  energy  transport  outside  the  plasma.  Energy  losses  include:  (1)  neutrons; 
(2)  charged  particles  that  do  not  deposit  their  energy  in  the  plasma.  Charged  particles  leaving  the  plasma  can 
be  heated  fuel  particles  or  energetic  fusion  products.  These  loss  particles  have  power  P|p;  and  (3)  unabsorbed 
radiation.  Bremsstrahlung  and  cyclotron  radiation  powers,  /V  and  P^ ,  respectively,  are  emitted  by  charged 
particles  due  to  interaction  with  magnetic  fields  present.  Thus  the  total  radiation  power  comes  at  the  expense 
of  a  reduction  in  kinetic  energy  of  the  particle  flow.  A  fraction,  7cy,  of  the  cyclotron  radiation  is  absorbed 
by  the  plasma,  while  bremsstrahlung  radiation  is  assumed  to  be  lost  from  the  plasma.  Although  the  injected 
power,  Pinj,  is  meant  to  heat  the  fusion  fuel  or  power  the  confinement  system,  it  might  also  contribute  to 
the  radiation  power  loss  as  a  side  effect.  For  example,  RF  heating  in  the  Field  Reversed  Configuration  would 
increase  particle  collisions  in  the  fusion  plasma,  increasing  the  bremsstrahlung  radiation. 
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The  required  injected  power,  Pi„j ,  is  found  by  performing  a  power  balance  on  the  fusion  plasma.  This 
power  can  be  imparted  by  charged  or  neutral  beam  injection,  Ohmic  heating,  or  radiowave  heating.  [7]  In 
order  to  maintain  a  constant  energy  in  the  plasma,  energy  must  be  injected  into  it  to  make  up  for  radiation 
and  particle  transport  losses.  Fortunately,  some  internal  heating  occurs  in  the  plasma  due  to  deposition  of 
kinetic  energy  of  charged  particles,  from  the  fusion  reaction.  If  this  internal  heating  is  sufficient  to  make 
up  for  the  losses,  the  plasma  is  said  to  be  ignited;  if  not,  external  heating  power  is  necessary.  Thus,  Pinj  is 
sum  all  the  power  loss  and  heat  addition  mechanisms  of  the  plasma.  The  heating  of  the  plasma  occurs  with 
efficiency,  ,  so  the  electrical  power  required  to  run  the  fusion  reactor  is  Pinj  / r/jn; . 

The  power  balance  for  the  plasma  is: 

Pf  +  Pinj  —  Pir  +  (1  —  1cy)P cy  +  Pip  +  (1  —  y)P t  (25) 


2.3  Propellant  Heating 


The  propellant  serves  as  a  coolant  that  absorbs  the  thermal  power  from  the  rocket  components.  A  minimum 
propellant  flow  is  set  that  adequately  cools  the  reactor,  mixing  chamber,  and  nozzle.  The  radiation  and 
neutron  power,  along  with  Ohmic  heating  ( Pohmie  =  I2R)  generated  by  the  magnetic  nozzle  (if  present)  and 
reactor,  and  the  heat  produced  by  the  injector,  constitute  the  thermal  power,  Ptherm- 

For  cooling,  the  propellant  must  flow  through  passages  near  the  hot  component  surfaces.  The  highest 
temperature  the  propellant  can  reach  is  limited  by  the  material  of  the  component  being  cooled.  To  maximize 
heat  transfer  and  total  energy  absorption,  the  propellant  should  enter  the  passages  at  the  lowest  temperature 
possible.  Assuming  steady  state,  the  thermal  power,  Ptherm,  absorbed  by  the  propellant  sets  the  necessary 


mass  flow,  according  to: 


Ptherm 
rnprop  —  .  . 

'tlMl  _  <» 0 


(26) 


where  hmat  and  h0  are  the  final  (materials  limited)  and  initial  enthalpies  of  the  propellant. 


2.4  Power  Generation  and  Demands 

A  turbine  extracts  work  adiabaticaliy  from  the  enthalpy  of  the  heated  working  fluid.  Conversion  of  this  work 
to  electrical  energy  occurs  at  an  efficiency  of  T)ten  ■  For  a  turbine-generator  system,  rjgtn  is  a  function  of  the 
temperature  difference  across  the  generator,  while  the  power  produced  is  a  function  of  the  enthalpy  difference 
across  the  turbine.  The  electrical  power  produced  must  meet  the  demands  of  the  payload,  propulsion  system 
components,  and  fusion  reactor.  The  heated  propellant,  y turboprop,  sent  to  the  turbines  is  sufficient  to  power 
the  electrical  generators  (0  <  ytUrb  <  1).  The  remainder  of  the  propellant,  (1  -  ytUrb)TTiprop,  bypasses  the 
turbines  and  is  sent  directly  into  the  mixing  chamber.  The  electrical  power  generated  is 

Ptl,gen  =  •7jeflT«urs(^ma*  ~  hturb)firlprop  (27) 

where  htUrb  is  the  enthalpy  of  the  propellant  exiting  the  turbines.  After  work  has  been  extracted  from  the 
heated  propellant  in  the  turbines,  the  propellant  is  mixed  with  the  escaping  plasma  particles  (with  power 

flp). 

The  power  required  for  a  magnetic  coil,  such  as  used  in  a  magnetic  nozzle  or  theta-pinch  coil,  depends 
on  the  required  magnetic  field,  the  resisitvity  of  the  material,  and  the  dimensions  of  the  magnet.  The  power 
necessary  to  produce  a  magnetic  field,  B,,  in  a  single-turn  solenoid  is  [7] 


( \fiog(a',0 ')) 


prfT\ 


where  pre,  is  the  resistivity  of  the  magnet,  A  is  the  volume  fraction  of  the  magnet  devoted  to  producing  the 
field  (around  90%),  is  the  inner  radius  of  the  coil,  and  g(a' ,0')  is  a  function  of  r1(  r2,  and  L, 


g(a',n  = 


2w(a'2  -  1) 


In  <*'  +  ■\A*/2  +  /?'2  j 

_  1  +  sqrtl  -4-  /?'J  J 


where  a1  and  are  defined  as  at'  =  r2/ri  and  (f  =  L/(2ri),  and  L  is  the  axial  length  of  the  magnet.  The 
outer  radius  r2  of  the  magnet  depends  on  the  radial  thickness  of  the  magnet,  which  in  turn  depends  on  the 
tensile  strength,  <rten ,  of  the  magnet  material. 


r2~r  i  = 


a**n  (  Trf  )  —  3 


Ohmic  heating  of  the  magnets  converts  the  electrical  power,  Pmagnet,  into  thermal  power  of  magnitude 
Pmagnet ■  This  contributes  to  the  total  thermal  power,  Pthtrm .  and  must  be  dissipated.  The  formula  for 
Pmagntt  can  be  used  to  find  the  power  demands  of  a  magnetic  nozzle  or,  in  the  case  of  a  FEC,  a  theta  pinch 
and  end  mirrors.  If  we  employ  superconducting  materials  for  our  magnets,  we  can  drop  the  resistivity  to  zero, 
and  thus  the  Ohmic  cooling  requirements  will  also  drop  to  zero.  Issues  that  remain  to  be  examined  include 
the  superconductive  temperature  threshold,  materials  properties,  and  effects  of  neutron  and  cosmic  radiation 
flux  on  superconductive  properties.  If  we  cannot  use  superconducting  materials,  another  alternative  is  to 
reinforce  the  magnets  with  high- tensile  strength  steel  bands,  in  effect  increasing  the  tensile  stength,  <rten, 
of  the  magnet.  However,  the  final  model  should  also  account  for  secondary  gamma  photons  produced  by 
neutron  flux  on  high-Z  materials,  such  as  steel  bands. 

The  payload  power  demands,  Ppavid,  are  dictated  by  the  payload  and  mission;  “housekeeping”  power  is 
around  40  KWe.[34j  The  total  quantity  of  required  electrical  power  is  taken  to  be  the  sum  of  the  magnet 
coils  and/or  capacitors,  magnetic  nozzle,  pumps,  power  injectors,  and  payload. 


,  =  ^  +  Preac  +  Pf 

Thni 


+  P. . 
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2.4.1  Reactor*  and  System-Ignition 

Ideally,  for  the  plasma  to  be  self-sustaining,  it  must  heat  itself  so  that  the  injected  power,  ,  is  minimized. 
To  accomplish  this,  we  must  maximize  the  fractions  of  charged  particle  and  cyclotron  power,  yPc  and  jeyPey 
respectively,  that  go  to  heating  the  plasma.  Defining  Qp<r  =  Pj  / Pin j  as  the  ratio  of  fusion  power  over 
injected  energy,  we  have  reactor  “breakeven”  when  Qp,r  =  l.[8] 

If  the  generated  electrical  power,  Pei, gen ,  is  insufficient  for  the  system  needs,  Pti,req,  an  external  source 
P,xt  =  Pti'ftq  -  Pei, gen  is  required.  When  the  system  requires  no  external  power,  we  can  establish  that  the 
system  is  “ignited”:  using  QPi,  =  Pat,g*n/Pext,  ignition  occurs  when  the  system  can  supply  its  own  electrical 
power  needs.  This  is  desirable  because  it  excludes  heavy  battery  or  solar  power  supply  systems  from  the 
rocket  mass,  improving  performance.  If  ignition  is  achieved  for  the  propulsion  system,  no  external  power  is 
required  (P«t  =  0),  so  Pei,gen  =  Pet, re,  and  Qp,t  -*  oo. 


2.5  Plasma  Mixing 

The  model  assumes  that,  in  the  “mixing  chamber” ,  the  fraction  of  heated  propellant  bypassing  the  turbine 
mixes  with  the  fraction  that  passed  through  the  turbine.  The  total  propellant  flow  is  then  heated  by  the 
energetic  charged  fusion  particles  to  a  uniform  temperature.  This  is  one  of  the  critical  assumptions  of  this 
model.  In  other  fusion  propulsion  studies,  propellant-plasma  mixing  has  been  identified  as  a  challenging 
area.  In  the  McDonell  Douglas-General  Atomics  study  [10],  the  chamber  in  which  plasmoids  mixed  with 
propellant  was  on  the  order  of  50  m  long  to  accommodate  adequate  thermal  transfer. 

The  propellant  mixes  adiabatically  and  uniformly  with  the  energetic  plasmarloss  particles  escaping  the 
fusion  plasma,  Pip.  Thus,  the  power  of  the  mixed  and  heated  propellant  becomes 

Pmix  =  Ptherm  +  Pip  ~  (32) 

Vgen 


Assuming  complete  ionization,  with  three  degrees  of  freedom  per  particle,  the  heat  capacity  of  the  particles 
can  be  found  by  Fermi-Dirac  statistics  [35,36]  to  be  Cp  =  5/2  R  J/mole-K,  where  R  =  8.3143  J/mole-K 
is  the  universal  gas  constant.  If  we  divide  both  sides  by  the  average  mass  of  a  mole  of  particles  (ions  or 
electrons),  we  have  Cp  in  units  of  [J/kg-K].  Assuming  steady  flow,  we  can  find  the  stagnation  temperature 
of  the  exaust  by  the  identity 


Pmix  — 


Cp 


(33) 


All  the  energy  of  the  propellant  is  random,  thermal  energy  instead  of  kinetic  energy  (stagnation).  The 
exhaust  is  then  directed  through  a  nozzle,  and  the  thermal  power  of  the  mixed  propellant,  Pmix,  is  converted 
to  jet  power. 


2.6  Thrust  Generation 

The  total  mass  flow  rate  exiting  the  rocket  is  the  sum  of  the  mass  flow  rates  of  the  fuel  and  the  propellant: 

rHou*  —  T7lprop  -f*  TTlfuel  (34) 

We  can  calculate  the  thrust  generated  by  the  propulsion  system.  Assuming  that  the  exhaust  temperature 
and  density  correspond  to  use  of  a  magnetic  nozzle,  we  can  use  Fig.  7  [31],  which  gives  us  state  ratios  for  a 
meridional  magnetic  nozzle.  We  need  the  ratios  of  inlet  temperature  to  throat  temperature  (Tm<x/T*),  and 
exit  velocity  to  throat  velocity  (v«*/v*).  In  Fig.  7  these  ratios  are  listed  as  (TJ=_4.o/TJ=o)  and  (vJ=4_Q/vz=0). 

By  Eq.  33  we  have  the  mixing  temperature,  Tm,x,  which  is  at  the  entrance  to  the  nozzle.  Dividing  Tm,x 
by  ( TmiX/T *)  gives  the  throat  temperature  of  the  exhaust.  At  the  throat,  the  exhaust  thermal  energy  is 
converted  to  kinetic  energy: 

fhoutCpT1*  —  —  rhovtv  (35) 
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Expansion  of  the  nozzle  after  the  choke  point  accelerates  the  particles  from  v *  to  vts.  Solving  for  i>*,  and 
multiplying  by  {vexh’)  from  Fig.  7,  gives  us  the  exit  velocity  of  the  exhaust.  Thrust  is  then  the  product  of 
the  exhaust  mass  flow  rate  times  the  exit  velocity. 

2.7  Component  Masses 

2.7.1  Coils  and  Magnetic  Nozzle 

The  reactor  mass  depends  entirely  on  the  type  of  reactor  employed:  if  an  FRC  is  used  for  fusion  power,  the 
theta  coil  and  mirror  magnets  will  be  the  dominant  reactor  mass;  if  a  DPF  is  used,  the  reactor  mass  should 
be  quite  low.  Similarly,  the  nozzle  mass  depends  on  whether  magnetic  coils  are  used  or  standard  nozzle 
materials. 

The  main  factor  in  determining  mass  of  magnetic  coils  is  the  tensile  strength  of  the  material  necessary 
to  withstand  the  forces  of  magnetic  pressure.  With  the  thickness  and  the  length  of  the  cylindrical  magnet, 
multiplying  with  the  magnet  material  density  gives  the  magnet  mass.  This  technique  will  work  for  the  theta 
pinch  and  mirrors  of  a  FRC,  and  the  magnetic  nozzle  of  the  FRC  and  DPF. 

In  Eqs.  (28),  (29),  and  (30),  the  radial  thickness  of  a  cylindrical  magnet  was  found  for  a  given  inner 
radius,  axial  length,  resistivity,  tensile  strength,  and  required  field  strength.  Given  the  density  of  the  magnetic 
material,  pmag,  the  mass  of  the  magnet  is  the  product  of  density  and  the  cylinder  volume: 

rrimagnet  =  *(2ri  Ar  +  (A  r)2)Lpmag  (36) 

where  Ar  =  r2  —  ri . 

2.7.2  Shielding 

The  mass  of  the  radiation  shielding  will  depend  on  the  duration  of  the  mission,  and  the  flux  and  type  of 
radiation  encountered.  It  is  assumed  that  the  reactor  is  the  only  source  of  man-made  radiation,  and  that 
natural  radiation  includes  trapped  charged  particles  (Van  Allen  belts),  Galactic  Cosmic  Radiation  (GCR), 
and  Solar  Particle  Events  (SPE’s  or  Solar  Flares).  While  the  positions  of  the  Van  Allen  belts  are  known, 
and  GCR  maxima  roughly  follow  an  eleven  year  cycle,  SPE’s  cannot  currently  be  predicted.  SPE  doses  can 
be  severe,  thus  it  is  prudent  to  include  on  board  an  SPE  ”  storm  shelter”.  For  a  quick  trip  to  Mars,  a  mass  of 
five  metric  tons  has  been  assigned  to  the  storm  shelter,  based  on  estimated  volume  requirements  and  shield 
densities[37]. 

Since  the  reactor  is  running  steady  state,  it  is  assumed  the  radiation  dose  rate  is  constant;  it  is  further 
assumed  that  the  total  dose  from  the  reactor  to  the  occupants  is  the  limiting  factor  in  determining  reactor 
shield  thickness.  Thus,  for  a  given  mission  duration,  the  shield  thickness  must  be  varied  to  prevent  the 
total  dose  from  exceeding  the  maximum  allowable  dose.  The  dose  depends  on  the  energy  of  the  attenuated 
neutrons  as  well  as  the  flux. 

Since  we  will  most  likely  be  employing  the  D(3He,p)a  reaction  in  our  reactor  [38],  there  will  be  an 
opportunity  for  parasitic  D-D(n)  reactions.  There  are  two  branches  to  D-D  [3]  with  almost  equal  cross 
sections,  as  shown  in  Table  2.  The  products  of  the  D-D  reactions  are  all  charged  particles  except  for  a 
2.45  MeV  neutron,  coming  from  the  D-D(n)  branch;  bremsstrahlung  radiation  produced  is  easily  shielded 
against. [39]  D-D  fusion  has  a  higher  <cv>  than  D-3He  at  almost  all  temperatures  [8],  so  the  D-D  source  of 
neutrons  is  not  likely  to  be  negligible.  Thus,  shielding  against  the  effects  of  neutrons  will  be  necessary. 

A  logical  choice  of  shielding  material  for  space-based  applications  is  lithium-hydride  (LiH)  [40];  it  is 
light-weight,  has  low-Z  material  to  thermalize  fast  neutrons,  and  is  a  well-  researched,  established  shielding 
material  whose  properties  are  known.  Radiation  shielding  is  a  complex  problem,  so  as  a  first  approximation 
I  have  assumed  a  thick  LiH  shield  (>  40  —  50  cm  for  »  1013  neutrons/cm2s)  that  attenuates  collided  neutrons 
completely  based  on  the  removal  cross  section,  and  those  that  pass  through  retain  all  their  energy1 .  Neutrons 

1  Using  this  thick  shield  approximation,  any  secondary  gammas  that  would  be  generated  would  exit  the  shield  at  the  same 
rate  as  the  uncollided  neutrons,  but  at  a  lower  energy  rate;  however,  since  we  are  using  a  low -Z  shielding  material,  secondary 
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entering  the  shield  collide  with  shield  material,  depositing  some  fraction  of  their  2.45  MeV,  perhaps  scattering, 
perhaps  being  absorbed.  By  assuming  straight  flux  attenuation  with  a  hydrogenous  thick  shield,  a  reasonable 
approximation  [41,42,43]  of  actual  shielding  results  is  possible  without  having  to  perform  numerical  multi¬ 
group  Monte  Carlo  calculations. 

The  source  of  the  neutrons  is  the  D-D(n)  reaction,  which  proceeds  at  a  reaction  rate  of  RR.D-D(n)- 
Assuming  isotropic  distribution  (the  neutrons  travel  in  all  directions  from  the  source  uniformly),  the  original 
flux  to  the  passengers  is 


We  will  use  a  spherical-shaped  shadow  shield  (known  as  a  “four-pi”  shield)  to  block  out  a  fraction,  Xcone . 
of  the  isotropic  neutrons,  taken  here  to  be  1/8  of  a  sphere.  The  shield  has  inner  and  outer  radii,  rj  and  rj, 
and  to  minimize  shield  mass  for  a  given  coverage,  we  wish  to  have  rj  as  close  to  the  source  as  possible.  As 
a  first  approximation,  the  flux  to  the  passengers  is  attenuated  by  the  thick,  hydrogenous  shielding  by  the 
relation 

$  =  $o  exp(— EArjh<e/<1)  (38) 

where  E  is  the  effective  cross  section  of  the  LiH  shielding,  given  by 


RRp-D(n)Vp 


™  'past 


£  =  nthieldivLi  +  ffH ) 

Pthield^AV ,  , 


MW  LiH 


( <*Li  +  &h) 


Solving  for  the  shield  thickness,  Ar,hieid,  we  can  calculate  the  shield  volume  by 


Vthield  —  Xcone  g  ^(•‘j  ri ) 


(39) 

(40) 


which,  when  multiplied  by  the  shield  density,  pun  =  0.82  g/cm3,  gives  the  shield  mass. 

We  establish  the  fined  neutron  flux,  <&,  and  therefore  the  shield  thickness  and  mass,  by  the  integrated  dose 
allowable  for  the  passengers/payload.  From  Chilton  [39],  ANSI  and  the  NCRP  listed  3.5  x  10~8  cSv-cm2 
for  2.5  MeV  neutrons  as  the  “prescribed  neutron  response  functional  datum  for  the  phantom-related  dose 
equivalent”;  thus,  given  the  total  time  of  exposure  to  the  flux,  $,  we  can  find  the  dose  to  the  passenger.  If 
we  limit  this  dose  to,  say,  5  cSv  annually3,  we  can  find  the  shield  thickness,  given  the  reactor  temperature 
and  thus  the  neutron  source  flux  <&o-  (See  also  40CFR61  for  exposure  limits  to  radiation  workers.) 


gammas  shouldn't  be  a  problem.  Design  considerations  for  teh  reactor  first  wall  and  surrounding  components  will  have  to 
include  secondary  gamma  production 

2  Specified  by  US  Nuclear  Regulatory  Commission  (NRC),  following  recommendations  from  the  National  Council  on  Radiation 
Protection  and  Measurements  (NCRP)  [37] 
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3  REACTOR  CONCEPTS 

3.1  Dense  Plasma  Focus 

The  Dense  Plasma  Focus  (DPF)  [44],  shown  in  Fig.  9  [45],  is  a  coaxial  electrode  configuration,  with  the 
anode  in  the  center  and  the  cathode  surrounding  it.  The  simple  geometry  and  compact  size  of  the  DPF 
make  it  attractive  for  fusion  space  propulsion  applications.  [46] 


Figure  9:  Dense  Plama  Focus  Schematic,  Mather  Configuration. 

The  model  used  for  the  DPF  is  similar  to  that  used  by  Leakeas.[47]  A  capacitor  bank  is  discharged  across 
the  electrodes,  ionizing  the  fusion  fuel  gas  between  them.  Current  flows  radially  between  the  electrodes, 
inducing  an  azimuthal  magnetic  field.  The  Jr  x  Bg  force  accelerates  the  ions  and  electrons  axially.  Motion 
of  these  particles  in  the  “rundown”  phase  rapidly  reaches  steady  state,  using  the  “snowplow”  model.[48] 

The  Dense  Plamsa  Focus  (DPF)  model  consists  of  two  parts:  rundown  and  pinch.  During  the  rundown 
phase,  the  gas  between  the  coaxial  Mather  [44]  geometry  electrodes  is  ionized  and  a  fraction  swept  forward 
(rundown  mass)  by  the  current  sheath.  The  electrode  lengths  are  set  such  that  when  the  current  sheath 
reaches  the  end,  the  current  delivered  by  the  capacitor  banks  is  at  a  maximum. 

In  the  pinch  phase,  the  trapped  plasma  is  constricted  by  the  magnetic  field  field  of  the  current  at  the  end 
of  the  electrodes  (maximum).  The  compressed  plasma  is  heated  like  an  ideal  gas  to  fusion  temperatures,  with 
a  Maxwellian  energy  distribution  assumed.  Fusion  reactions  occur  at  a  rate  corresponding  to  the  density  and 
temperature  of  the  fuel  plasma,  with  the  accompanying  energy  release  going  to  heat  the  rocket  propellant. 

Theory  behind  the  DPF  model  is  presented  here.  Appendix  B  gives  some  sample  program  input/output 
and  the  source  code  to  the  DPF  circuit-equivalent  program,  cir.c. 


3.1.1  Plasma  Rundown 

During  the  rundown  phase,  the  DPF  is  approximated  by  an  LC  circuit,  depicted  in  Fig.  10. 

The  time  rate  of  change  of  the  inductance,  L,  is  proportional  to  the  plasma  rundown  velocity,  experimen¬ 
tally  observed  to  be  roughly  constant.[49]  For  devices  with  bank  energy  W  «  100  kJ,  the  rate  of  inductance 
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Figure  10:  DPF  Rundown  Circuit 


change  L  20  mO.  (H/s=fi.)  For  a  coaxial  transmission  line,  the  inductance  change  is 

(4I) 

where  vrun  is  the  rundown  velocity,  and  re  and  ra  are  the  cathode  and  anode  radii,  respectively. 

The  circuit  equation  is 

U  -  LI  -  LI  =  0  (42) 

L  =  Lq  +  Lt  (43) 

where  U  is  the  capacitor  bank  voltage  and  I  is  the  current  delivered.  Pinch  and  line  resistance  are  neglected 
to  simplfiy  the  circuit.  If  we  hold  the  capacitor  bank  energy,  Wt  the  initial  rate  of  current  change,  Ip,  and 
L  constant,  we  can  solve  this  circuit  in  dimensionless  form: 


where 


u'-r  - (if0*  +  nT  =  o 

h 


( wi0L*y /» 

r  =  — — 

(' WIo/L )1/3 

r  =  — U- 

{W/I3L)V  3 

In  terms  of  U * ,  the  circuit  equation  reduces  to 

U*  +  lff0'  +  lff{U5+n0*  =  0 

which  can  be  solved  with  a  series  solution  of  the  form 


f7*(f)  = 

n=0 


(44) 

(45) 

(46) 

(47) 

(48) 

(49) 


a0  =  US 

ai  =  0 

_ _ QnUp _ ar»+l(u  4~  1) 

an+2  “  2(n  +  l)(n  +  2)  tf0*(n  +  2) 

This  solution  is  a  polynomial  in  t *;  the  recursion  relation  for  coefficients  A?  and  above  is  also  given. 


(50) 

(51) 

(52) 
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The  voltage  on  the  capacitors  is  proportional  to  the  charge,  and  the  current  delivered  as  a  function  of 
time  is  simply  the  negative  time  derivative  of  the  capacitor  charge. 

r  =  --jLir  (53) 

Shown  in  Fig.  11  are  the  dimensionless  voltage  and  current  profiles  as  a  function  of  time  during  the 
rundown,  which  follow  expected  trends.  If  we  choose  Uq  =  2.12  [49],  we  get  the  maximum,  optimum  current 
of  I*  =  0.64  with  a  rise  time  of  t*  =  1.45. 


Figure  11:  DPF  Dimensionless  Voltage  and  Current  Profiles 


3.1.2  Pinch 


The  pinch  is  approximated  by  a  cylindrical  plasma,  whose  dimensions  remain  constant  during  the  stable 
pinch  lifetime.  When  the  pinch  is  in  equilibrium,  the  magnetic  pressure  at  the  plasma  edge  must  balance 
the  kinetic  pressure  of  the  plasma  (Bennett  relation), 


8 


=  nkT 


(54) 


where  rp  is  the  pinch  radius,  n  =  n<  +  ne  is  the  particle  number  density  in  the  pinch,  and  kT  is  the  particle 
energy.  The  pinch  becomes  unstable  when  current  is  no  longer  delivered  from  the  capacitors  (U  =  0),  and 
thus  the  magnetic  field  drops  to  zero. 


Pinch  Dimensions.  The  pinch  radius  during  the  pinch,  rp,  is  taken  to  be  [1] 

rp  =  r0/kr  (55) 

where  r0  is  the  original  pinch  radius  and  kr  is  the  compression  ratio.  A  compression  ratio  of  kr  =  10  is 
conservative.  We  take  the  initial  pinch  radius  to  be  the  anode  radius;  it  is  actually  less  than  the  anode 
radius,  but  this  is  conservative  for  a  given  compression  ratio. 


r0  »  r„  (56) 

While  compression  ratios  of  kr  =  10  are  common  for  Z-pinches,  DPFs  can  have  kr  =  100. [50] 

From  experimental  scalings,  we  take  the  pinch  length  to  be  approximately  the  difference  between  the 
cathode  and  anode  radii. 

IpBsr e-  ra  (57) 
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Pinch  Lifetime.  To  simplify  the  DPF  pinch  physics,  we  assume  that  the  pinch  is  stable  for  the  time  it 
takes  the  capacitors  to  discharge  the  rest  of  their  energy  at  the  maximum  current  and  voltage  at  maximum 
current  (dc).  Theoretically,  this  estimate  of  pinch  lifetime  is  a  best-case  scenario.  The  maximum  current 
flows  through  the  pinch,  and  the  pinch  is  stable  for  as  long  as  the  capacitor  bank  has  energy  stored.  Thus 
any  subsequent  calculations  of  energy  yield  from  this  pinch  will  represent  an  upper  limit.  However,  compared 
to  empircal  scalings  to  high  bank  energies,  the  dc  assumption  gives  conservative  pinch  lifetime  predictions. 

The  capacitor  bank  energy  devoted  to  rundown  of  the  plasma  is 


Wrd  =  U(t)I(t)dt 
Jo 


where  tmaa  is  the  time  at  which  the  current  reaches  a  maximum.  The  remaining  bank  energy  is  devoted  to 
the  pinch: 

Wpi„  -  W  -  Wrd  (59) 


The  stable  pinch  lifetime  is  then  taken  to  be 


U(t  max  )i(t  max) 


Pinch  Temperature.  When  we  consider  fusion  heating  and  radiation  cooling  of  the  pinch  plasma,  we  no 
longer  have  a  pinch  in  equilibrium.  For  simplicity,  it  is  assumed  that  the  pinch  dimensions  remain  constant 
over  the  stable  pinch  lifetime.3  The  initial  plasma  temperature  is  determined  by  the  Bennett  relation  (Eq. 
54),  and  the  subsequent  plasma  temperature  is  determined  by 

—  NkT  =  U  {tm  ax')  I  it  max')  +  Pjution  Phr  (61) 

where  N  is  the  particle  inventory  in  the  pinch,  and  Pjution  and  /V  are  the  fusion  and  Bremsstrahlung 
powers  that  heat  and  cool  the  plasma. 


3.1.3  Bremmstrahlung  Production 
The  Bremsstrahlung  emission  spectrum  is  [1] 


12.40  \ 

V  Akre  ) 


where  A  (A)  is  the  photon  wavelength.  The  “Gaunt  factor”  g  approaches  2^/Z/r  at  high  plasma  temperature. 
The  peak  Bremsstrahlung  emission  wavelength  is 

Ama*  =  (63) 


The  total  Bremsstrahlung  power,  Py.  (W/m3),  is  then 

*  = 

=  5.35  x  10 -^Z'jjnKkT)3  (64) 

3  In  fact,  the  pinch  radius  should  not  remain  constant  if  there  is  net  heating/cooling  of  the  plasma,  but  the  only  alternative 
to  this  simplifying  assumption  is  a  2-D  time-dependent  MHD  code. 
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where  kT  (keV)  is  the  plasma  Maxwellian  temperature,  and  nt  (m-3)  is  the  electron  number  density.The 
effective  atomic  number,  Zaj{,  is  defined 


Z'!I 


np  +  ri3jfe22 
n, 

5 

2 


(65) 


where  the  deuterium  density,  np,  is  equal  to  the  helium-3  density,  n3#e.  Nitrogen  plasma,  which  has  a 
higher  Z  so  gives  a  larger  Bremsstrahlung  contribution,  but  is  not  a  fusion  fuel,  but  has  Z*U  =  7. 

We  approximate  the  Bremsstrahlung  yield,  E^r  (J)  by 


Et  r  =  Vp 


(66) 


where  Vp  is  plasma  volume  (m-3),  rp<n  (s)  is  stable  pinch  time,  lPi„  (m)  is  pinch  length. 


3.1.4  Issues 

There  are  issues  to  be  addressed  and  assumptions  to  be  validated  before  the  DPP  can  be  employed  in  a 
propulsion  concept. 

•  Particles  in  the  pinch  are  thermalized,  their  energy  distribution  is  Maxwellian,  and  that  Maxwellian 
reaction  rate  parameters  (<<tv>)  may  be  used.  In  fact,  the  cross  sections  in  the  pinch  may  not  be 
entirely  (or  at  all)  Maxwellian,  but  some  combination  of  beam-target  and  Maxwellian.  Until  the  true 
temperature-cross  section  relationship  in  the  pinch  is  known,  there  will  be  a  degree  of  uncertainty  in 
specifying  the  fusion  output  power. 

•  The  kT  oc  /J  scaling  relation  is  valid  up  to  fusion  temperatures.  More  accurate  models  of  the  Dense 
Plasma  Focus  exist  [51,52,53],  including  time-dependent,  three-fluid  MHD  computer  simulations.  This 
model  of  the  DPF  used  is  rudimentary,  but  scaling  of  the  plasma  focus  device  has  not  yet  been 
experimentally  verified  for  currents  greater  than  /  «  1  MA  or  energies  greater  than  kT  ss  1  keV.[54] 

•  Pinch  lifetimes,  which  are  on  the  order  of  150  ns  [55],  can  be  made  orders  of  magnitude  longer,  through 
the  use  of  embedded  magnetic  fields  that  prevent  instability  modes.  Alternately,  capacitors  or  other 
power  supplies  can  be  made  to  fire  at  very  high  rates  of  fire  ( v  «  104  s-1)  to  compensate. 

•  The  pinch  is  a  cylindrical  plasma  with  no  internal  magnetic  field,  whose  initial  dimensions  are  measured; 
the  nature  of  the  time  rates  of  change  of  the  pinch  energy  and  dimensions  is  calculable.  The  short 
time  scales  and  complicated  plasma  behavior  of  the  pinch  make  accurate  theoretical  or  experimental 
models  difficult  to  realize. 

•  Fusion  fuel  and  charged  fusion  products  are  confined  during  pinch  lifetime;  it  follows  that  all  fusion 
energy  to  charged  particles  is  confined  (absorbed)  in  the  pinch  plasma  during  pinch  lifetime.  Experi¬ 
mental  observations  indicate  that  large  particle  emissions  occur  during  plasma  focus  pulses  [56,57,58], 
but  this  model  assumes  that  particle  transport  has  been  minimized  for  each  pulse. 

•  Electrode  materials  can  be  made  to  withstand  the  pulsed,  multi-mega  ampere  currents  delivered. 

•  Power  supply  specific  masses  (energy  per  unit  mass)  can  be  enhanced  so  supply  masses  are  not  pro¬ 
hibitively  large  for  space  applications. 
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Figure  12:  Field-Reversed  Configuration 


3.2  Field-Reversed  Configuration 

Fig.  12  shows  a  cross  section  of  a  Field  Reversed  Configuration. [59]  It  is  a  member  of  the  Compact  Toroid 
(CT)  family,  characterized  by  torus-shaped  plasmas  without  enclosing  magnetic  coils  or  core.  The  FRC  is 
stretched  along  the  center  axis  (prolate),  has  little  or  no  toroidal  field,  and  has  high  beta,  0.  The  FRC  has 
many  appealing  features  for  space  propulsion: 

•  The  design  of  the  FRC  burn  reactor  is  a  simple  cylindrical  tube,  simplifying  construction  and  mainte¬ 
nance,  enhancing  reliability,  encouraging  modular  design,  and  reducing  cost. 

•  The  FRC  reactor  does  not  have  a  solid  core,  as  does  the  tokamak,  and  translation  of  the  FRCs  has 
been  demonstrated.[60]  Thus  it  is  possible  to  separate  the  formation  and  burn  chambers  to  optimize 
the  characteristics  of  both. 

•  The  FRC’s  open  field  lines  act  as  a  natural  diverter  to  exhaust  burnt  fuel  (ash)  and  heat  the  propellant. 
The  open  field  lines  that  surround  the  FRC  would  permit,  in  steady  state  operation,  the  flow  of 
propellant  to  mix  with  the  hot  plasma  boundary  and  heat  the  propellant. 

•  Since  and  the  open  field  lines  travel  out  the  end  of  the  theta  pinch  coil,  we  can  take  advantage  of  the 
simple  geometry  of  the  FRC  by  having  several  FRC  burn  chambers  in  line,  with  the  open  field  lines 
flowing  through  each  of  them. 

•  FRC’s  seem  ideally  suited  for  rocket  use  in  conjunction  with  a  magnetic  nozzle.  The  field  free  core  and 
sharp  field  reversal  more  easily  permit  charged  particles  to  detach  from  the  magnetic  field  lines  and 
escape  from  the  nozzle. [7, 59] 

Although  FRCs  demonstrate  macrostability,  they  are  subject  to  microinstabilities  and  have  life-times 
on  the  order  of  300  /js.[61]  The  n  =  2  rotational  instability  has  recently  been  eliminated  by  using  weak 
multipole  fields  after  FRC  formation.[59]  Ultimately,  it  is  hoped  that  the  FRC  will  burn  in  steady-state;  the 
FORTRAN  code  assumes  steady  state  FRC  operation  and  refueling,  but  a  more  conservative  case  would  be 
to  find  the  pulse-averaged  fusion  power,  as  is  done  for  the  DPF. 

FRC  Formation.  FRCs  are  commonly  generated  by  field-reversal  of  a  single-turn,  elongated  solenoid 
(theta  pinch).  FRC  generation  can  be  described  in  four  stages  [61]:  pre-ionization,  field-line  connection, 
heating,  and  containment. 

Methods  for  preionization  include  Z-pinch  injection  and  conical  theta-pinch  discharge  at  the  ends,  but 
the  most  common  method  is  discharging  a  theta  pinch  coil.  The  fill  gas  in  the  formation  chamber  is  turned 
into  a  cold  plasma  during  the  preionization  stage,  freezing  the  bias  magnetic  field  direction.  Once  the  bias 
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is  established  in  the  plasma,  the  current  in  the  theta  coil  is  reversed,  reversing  the  direction  of  the  magnetic 
field  near  the  coil  surface.  The  field  direction  in  the  center  remains  the  same  as  the  original  bias.  The 
magnetic  field  lines  are  connected  at  the  ends  of  the  theta  pinch  coil,  containing  the  fusion  plasma.  Field 
line  reconnection  is  traditionally  induced  by  passive  or  independent  fast-trigger  end  mirrors.  Finally,  the 
closed  field  lines  contract  axially  and  come  to  equilibrium  with  the  fields  from  the  end  mirrors;  adiabatic 
contraction  heats  the  plasma  to  the  burn  temperature. 


Reactor  Scaling.  Studies  concerning  FRC  lifetime  scaling  [62]-[80],  reflect  the  desire  for  steady  state 
operation.  Particle  transport  out  of  the  closed-field  line  region  is  a  prominent  factor  in  limited  configuration 
lifetimes.  A  typical  particle  is  contained  within  the  separatrix  for  an  average  time,  tjv  ,  which  is  defined  as 
the  particle  inventory  divided  by  the  loss  rate.  It  may  not  be  possible  to  contain  particles  absolutely,  but  it 
is  desirable  to  make  r n  as  large  as  possible.  As  we  increase  reactor  size,  temperature,  and  particle  density 
to  reactor  conditions,  an  important  issue  is  to  see  how  the  confinement  time  scales  with  these  parameters. 

A  critical  issue  to  the  validity  of  the  FRC  thruster  model  is  scaling  to  higher  temperatures.  Empirical 
scaling  laws  have  been  fitted  to  experimental  data  [65,66,71,72,76,78,79,80],  but  there  is  indication[59]  that 
the  results  of  particle  and  flux  confinement  and  configuration  lifetime  are  very  dependant  on  seemingly 
unrelated  factors  as  fill  gas  pressure,  pre- ionization  technique,  and  reactor  dimensions. 

At  this  stage  of  FRC  research,  effort  is  directed  towards  demonstration  of  the  concept  through  prolonging 
confinement  and  lifetime. [59]  Ultimately  the  FRC  reactor  would  operate  at  steady  state  conditions,  with  fuel 
injection  to  compensate  for  particle  loss  and  fuel  burn-up. [38, 81]  There  are  currently  no  experimental  data 
of  FRC  operation  at  reactor-like  conditions,  and  so  the  model  relies  upon  theoretical[67,68,69,70,73]  and 
experimental[77]  predictions. 

Steinhauer[77]  suggests  that  at  reactor  conditions  particle  confinement  would  scale  classically  with  tem¬ 
perature  and  linearly  with  size  (effective  diffusivity  inversely  proportional  to  size).  Krall’s[69]  proposed 
scaling  law  is  based  on  low-frequency  drift  wave  transport,  agreeing  with  other  theoretical  and  experimental 
results  in  that  as  the  collision  frequency  decreases,  energy  and  particle  loss  also  decreases. 

Miley’s  Velocity  Space  Loss  Sphere  (VSLS)  model[67,68]  represents  the  upper  bound  of  particle  con¬ 
finement  in  the  FRC.  It  assumes  that  the  dominant  loss  mechanism  is  due  to  scattering  into  a  region  of 
( H,P» )  space,  or  velocity  space,  where  particles  are  not  contained;  other  mechanisms,  such  as  instability  or 
cross-field  diffusion  are  negligible,  i.e.  have  been  solved.  ( H  and  Pg  are  the  total  energy  and  canonical  angu¬ 
lar  momentum  of  the  particle  considered,  respectively.)  VSLS  confinement  represents  the  best  confinement 
possible  because  it  implicitly  assumes  having  overcome  current  engineering  challenges.  A  fraction,  between 
one  and  ten  percent,  of  the  VSLS  scaling  law  is  used  to  represent  particle  confinement.  This  allows  the 
model  to  provide  optimistic  forecasts  for  the  FRC’s  potential,  but  prevents  the  model  from  overestimating 
the  the  FRC’s  ability. 

Table  3  summarizes  features  of  these  four  models. 


Table  3:  FRC  Particle  Confinement  Time  Scaling:  rjy  ~  If  xb,r‘ndTe 


Scaling 

a 

b 

c 

d 

e 

TRX-2  Experiment [59 ,77] 

0 

1.0 

1.0 

0.9 

1.5 

VSLS  Theory[73] 

0 

2.7 

1.8 

-1.0 

0.2 

Krall  Theory[69,70] 

0 

2.0 

2.0 

0.5 

-0.5 

Lower  Hybrid  Drift  Theory  [59] 

0.6 

3.3 

2 

0.5 

-0.7 
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3.2.1  Power  Balance 


The  FRC  is  assumed  to  run  in  steady  state,  so  the  temperature  in  the  plasma  remains  constant.  The  power 
balance  for  the  plasma  may  be  written 

~^Epla$  —  Pinj  +  Pf  ~  Pn  ~  Pip  —  (1  —  1cy)Pcy  ~  Pir  ~  0  (67) 

where  yCy  is  the  fraction  of  synchrotron  radiation  power  absorbed  in  the  plasma.  Assuming  Maxwellian  cross 
sections  apply,  the  fusion  power,  Pj ,  is  defined  as 

Pf  —  y  ]  RRj  Efj  Vpiat  (68) 


where  the  fusion  power  is  divided  up  between  neutrons,  Pn,  and  charged  particles,  Pe.  For  D-3He,  PJPj  as 
0.97. [8]  The  the  plasma  volume,  Vpiat,  is  assumed  to  be  defined  by  the  separatrix  boundries,  Vpia,  =  itr2tl ,, 
where  r,  and  l,  are  the  separatrix  radius  and  length,  respectively. 

For  a  FRC,  the  plasma  beta  is[59] 


<  0  > 

x. 


r./rc 


(69) 

(70) 


This  can  be  used  to  calculate  the  bremsstrahlung  and  synchrotron  radiation,  /V  +  Pcy,  described  in  Eq.  (4). 
The  particle  transport  power  out  of  the  FRC  plasma,  Pip,  can  be  expressed 


P,p  = 


NjkTjjl  +  f,rt 

Tft 


+  (1  -  y)Pc 


(71) 


where  y  is  the  fraction  of  the  fusion  charged  particle  power,  Pe,  that  is  absorbed  in  the  plasma,  and  rN  is 
the  particle  confinement  time,  scaling  for  which  is  found  in  Table  3. 


3.2.2  Issues 

Restart  Capability.  For  a  deep  space  mission,  it  may  be  necessary  to  coast  for  long  periods  of  time  or 
power  down  the  reactors.  Given  the  current  difficulty  in  achieving  ignition  in  D-T  fusion  reactors,  as  well  as 
the  higher  plasma  temperatures  necessary  for  D-3He  fusion,  the  issue  of  remote  restart  capability  deserves 
consideration.  Several  options  exist: 

1.  If  employing  multiple  fusion  reactors  to  produce  thrust,  one  reactor  could  be  kept  “lighted”  at  all  times 
to  act  as  a  “pilot” [82]  for  the  other  reactors.  Fusion  fuel  consumption  might  prohibit  keeping  even  one 
ignited  fusion  reactor  burning  at  all  times. 

2.  Cold  restarts  could  be  made  at  all  times  with  use  a  conventional  space-based  nuclear  fission  reactor, 
such  as  a  SP- 100(83],  as  a  starter  motor.  If  this  approach  is  taken,  additional  shielding  requirements 
would  have  to  be  taken  into  account. 

3.  Missions  within  the  solar  system  might  permit  use  of  solar  panels  to  trickle  charge  high  specific-energy 
capacitor  banks;  these  would  be  discharged  to  ignite  a  chemical[84]  or  D-T  fusion  “ramp”  [85],  which 
would  in  turn  ignite  the  D-3He  reactor(s).  Mass  of  necessary  solar  panels  and  capacitor  banks,  if 
necessary,  might  be  prohibitive. 

Quenching.  If  the  open  field  lines  are  to  be  used  to  carry  the  propellant  around  the  FRC,  where  mixing 
with  the  plasma  edge  heats  the  propellant,  too  much  propellant  might  overcool  the  FRC,  initiating  instabil¬ 
ities  or  taking  the  FRC  out  of  the  ignited  regime.  Thus,  there  might  exist  a  maximum  amount  of  propellant 
allowable  to  be  heated  directly  by  the  FRC,  depending  on  the  FRC  dimensions  and  stabilization  methods. 
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Unburnt  Fuel  Recovery.  Due  to  the  scarcity  of  3He  sources,  and  the  high  particle-loss  rate  of  an  FRC, 
unburnt  3 He  should  be  recovered  before  it  is  exhausted  from  the  nozzle;  recovery  apparatus  will  add  to  the 
system  mass,  but  it  might  be  less  expensive  to  the  mission  than  exhausting  costly  fuel. 

One  possible  method  to  retain  unburnt  fuel  would  be  to  convert  the  entire  system  to  nuclear-electric 
propulsion  instead  of  nuclear-thermal:  using  the  Venetian  blind  type  direct  conversion  (dc)  method[8],  fusion 
power  could  be  converted  to  electricity,  which  would  in  turn  power  a  magneto-plasma-dynamic  (MPD) 
or  arcjet  propulsion  device.  Recovery  of  the  3He  from  the  dc  device  would  then  be  simplified.  However, 
the  multiple  conversions  required  for  nuclear-electric  propulsion  would  likely  be  less  efficient  than  direct 
nuclear-thermal. 

The  objective  of  recovering  the  unburnt  fuel  is  to  reinject  it  into  the  plasma  for  reburning,  to  increase 
the  efficiency  with  which  the  3He  is  used.  If  the  reactor  type  is  an  FRC,  multiple  steady-state  CTs  could  be 
maintained,  separated  by  mirror  coils,  analogous  to  the  series  tandem-mirrors  in  the  SAFFIRE  concept[86]; 
the  open  held  lines  would  flow  through  each  chamber  in  series,  carrying  with  them  unburnt  fuel  and  ash  (burnt 
fuel);  these  open  field  lines  act  as  natural  diverters [59],  and  in  successive  chambers  the  particle  flow  could  be 
skimmed  off,  the  3He  separated  and  reinjected.  Using  this  method,  propellant  would  only  flow  around  the 
last  FRC,  as  (1)  no  unburnt  fuel  diverting  would  be  possible,  and  (2)  propellant  injected  earlier  in  the  series 
would  pollute  the  separation  process.  Since  translation  of  FRCs  has  already  been  demonstrated[60,79],  one 
formation  chamber  could  be  used  for  all  FRCs,  which  would  then  be  sent  to  their  respective  burn  chambers. 
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4  Example  Problem 


An  example  problem  illustrates  the  main  features  of  this  model,  using  the  DPF  as  the  fusion  reactor. 

We  start  with  a  capacitor  bank  energy  W  =  10  MJ,  an  initial  rate  of  current  change  Io  =  1  MA/^s,  and 
a  rate  of  inductance  change  L  =  20  mfi.  (1  H/s  =  1  Q)  This  gives  the  characteristic  voltage,  Au,  current, 
A/,  and  time,  At: 

Au  =  (Wi0L2)l/3  =  159  kV  (72) 

At  =  (Wi0/L)1/3  =  7.94  MA  (73) 

At  =  (W/%L)1/3  =  7.94 (74) 

Using  /*  =  {/’=  0.64  and  t*  =  1.43,  the  maximum  current,  voltage  at  maximum  current,  and  the  rise  time 
are 

Im  =  5.08  MA  (75) 

Um  s  102  kV  (76) 

tr  =  11.4  /is  (77) 


The  energy  of  the  capacitor  banks  is  divided  between  the  rundown  and  the  pinch:  W  =  Wrj  +  WPi„ 
The  fraction  of  bank  energy  devoted  to  rundown  is 


Wrd  f1A3 

W  J0 


U'Vdf  =  0.89 


(78) 


so  the  rest  goes  to  the  pinch: 


Wpin/W  =  0.11 


(79) 


This  compares  favorably  with  experimental  results,  e.g.  Gates[87]  gets  Wpin/W  ss  0.12  —  0.24. 

We  assume  the  pinch  is  stable  for  the  time  it  takes  the  capacitor  banks  to  deliver  the  remainder  of  energy 
at  voltage  U  —  Um  and  current  I  =  Im.  Although  this  might  not  be  accurate  (e.g.  pinch  might  be  stable 
only  as  long  as  current  increases),  the  assumption  conserves  capacitor  bank  energy,  and  gives  a  conservative 
pinch  lifetime.  Noting  that 

AuAiAt  =  W  (80) 


the  pinch  lifetime  is 


Wpin 

UmIm 

0.11  W 

^0.64At/)(0.64A/) 


(0.64)2 
2.13  na 


(81) 

(82) 

(83) 

(84) 


This  order  of  magnitude  compares  with  Herold[88],  who  gets  rpin  =  11  ps  at  /  =  20  MA. 
The  initial  plasma  temperature  is  set  by 


where 


uiT  /io/2 

(85) 

n  =  N/irr2inlpin 

(86) 

rpin  =  ro/kr 

(87) 

Ipin  re  ~  ra 

(88) 
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Taking 


re  =  8  cm  (89) 

r„  =  3  cm  (90) 

N  =  6  x  1019  (91) 

I  =  Im  (92) 

the  plasma  radius,  volume,  density,  and  initial  temperature  are 

rpin  =  0.428  mm  (93) 

Vpi„  =  2.88  x  10"2  cm3  (94) 

n  =  2.08  x  1021  cm"3  (95) 

kT0  =  6.71  keV  (96) 


For  now  we  assume  that  the  fraction  of  fusion  fuel  burnup  during  the  pinch  lifetime  is  small  compared 
to  the  fuel  inventory,  so  that  we  can  treat  the  fuel  density  as  constant.  (We  check  this  after  calculating 
the  fusion  yield.)  Rather  than  doing  severed  iterations  by  computer,  we  get  a  rough  idea  of  the  fusion  and 
radiation  yield  by  doing  one  Runge  Kutta  (RK)  iteration  by  hand. 

Whereas  the  Euler  method  integrates  the  temperature  by 


kT(ta  +  h)  »  kT(to)  +  gh 

(97) 

where  g=dkT/dt  is  the  local 

“slope” ,  RK  integrates  using  the  weighted  average  of  four  Euler  slopes: 

*  =  ikT 

kT0 

(98) 

*  =  ikT 

kTo+gih/2 

(99) 

*  ■  itT 

^To+^a^/2 

(100) 

«  •  itT 

(101) 

The  RK  estimate  is  then 

kT(t0  +  h)  *  kT(t0)  +  {j  +  j  +  j  +  j)h 

(102) 

We  start  with  the  the  time  rate  of  change  of  the  plasma  energy: 

±kT  =  ( Umlm  +  Pj  -  Pir) 

(103) 

The  power  delivered  from  the  capcitor  banks  is 

UmIm  =  (102  kV)(5.08  MA) 

(104) 

=  0.516  MJ//1S 

(105) 

The  fusion  heating  power  is 

Pj  = 

WjVpi^  <<TV> 

(106) 

ss 

(18.3  Me V)(  1.602  x  10" 

13  J/MeV)(2.88  x  10~2  cm3) 

(107) 

(2.08  x  1021  cm"3)2 .  3  .  . 

^ - - - —{<  trv  >  cm3/s) 

4 

= 

9.14  x  1016  <  <rv  >  MJ/jts 

(108) 
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Similarly,  the  radiation  cooling  is 


Pbr  =  AZ'jfV^n'VkT  (109) 

=  (5  x  10_31)(2.5)(2.88  x  10"2  cm3)(2.08  x  1021  cm“3)2(fcT  (keV))1/2  (110) 

=  0.156/bT1/2  MJ//is  (111) 

Our  function  reduces  to 

JtkT  =  6~x~10^(Q-516  +  9'H  X  1Ql6  <<TV>  -O-1^2)  (112) 

keV  106  J 
1.602  x  10“16  J  MJ 

=  (64.4  +  9.51  x  1018  <  (tv  >  -16.2ibT1/2)  keV/ps  (113) 


Table  4  lists  some  D-3He  <  crv  >  values  for  convenient  plasma  temperatures. 


Table  4:  Useful  <  crv  >  values 


kT  (keV) 

<  crv  >  (cm3/s) 

1 

2.1  x  10- 

10 

3.0  x  10"19 

100 

1.3  x  10“16 

1000 

2.5  x  10" 16 

5000 

2.0  x  10“16 

Using  h  =  rpin  =  2.13  /is,  we  approximate  the  temperature  change  as  shown  in  Table  5 


Table  5:  RK  Temperature  Change  Approximation. 


t 

kT  (keV) 

g  =  dkT/dt  (keV/^s) 

to 

6.71 

12.0 

=  l/i 

to  +  h/2 

19.5 

21.1 

=  92 

to  +  h/2 

29.1 

96.0 

=  93 

to  +  h 

211 

2040 

=  94 

At  the  end  of  the  pinch  lifetime  (f  =  0  +  rpin),  the  plasma  temperature  is  approximately 

kT(rPin)  ss  6.71(keV)  +  ^(12.0  +  2(21.1)  +  2(96.0)  +  2040)(keV//is)(2.13  /is)  (114) 

=  818  keV  (115) 

We  take  the  average  temperature  to  be  the  temperature  at  t  =  Tpin/ 2,  which  is  kT  =  412  keV.  We  then 
approximate  the  fusion  and  bremsstrahlung  energy  yields  by  multiplying  the  pinch  lifetime  by  the  powers  at 
the  average  temperature: 


Ej  Tpin  P f(kTav )  (HO) 

Ehr  ss  rpi„/V(*Ta„)  (117) 

This  gives  Ej  «  53.1  MJ,  and  Eb r  w  6.74  MJ. 

We  now  can  check  our  assumption  that  the  fusion  fuel  burned  is  much  less  than  the  original  inventory. 
The  number  of  fuel  ions  burned,  IV/,  is  twice  the  pulse  fusion  energy  divided  by  the  energy  per  fusion 
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reaction: 


w  - 

2(53.1  MJ) 

"  (18.3  MeV)(1.602  x  10"13  J/MeV)(6  x  1019) 

=  0.603 


(118) 

(119) 

(120) 


So  we  have  Nj  <  N,  as  we  should,  but  we  don’t  have  N/  <<  N,  which  we  need  for  our  assumption  that 
n  as  constant  during  the  pinch  lifetime.  However,  once  the  fuel  density  drops  enough  to  affect  the  reaction 
rate,  most  of  the  fuel  has  burned  up  already.  Thus  for  this  example  problem,  we  are  relatively  safe  as  long 
as  Nj  <  N. 

If  we  fire  at  a  repetition  rate  of  v  —  1  s- 1 ,  our  time  averaged  fusion  and  radiations  powers  are 


Pj  =  Eju  (121) 

P'ir  =  Eb rU  (122) 


We  set  the  propellant  flow  rate  to  absorb  thermal  power  produced  in  the  rocket.  The  thermal  power  is 
the  sum  of  the  input  (capacitor)  and  output  (fusion,  radiation,  nozzle  ohmic  heating)  powers.  Here  we  take 
this  to  be 

PtH=(W  +  Ef+Eir)u+Pn„  (123) 

Pno,  depends  on  the  nozzle  cooling  requirements,  i.e.  the  current  flowing  through  the  magnetic  nozzle,  which 
depends  on  the  magnetic  field  required  for  the  degree  of  ionization  of  the  propellant.  For  now  we  assume 
that  Pnoi  is  negligible  compared  to  the  radiation  and  injected  powers,  and  will  check  this  later. 


Pth  «  (W  +  E}+Ebr)u  (124) 

=  (10  +  53.1  +  6.74)  MW  (125) 

=  69.8  MW  (126) 


The  propellant  flow  required  to  absorb  this  thermal  power  depends  on  the  temperature  increase  we  can 
give  the  propellant,  which  depends  on  the  temperature  limitations  of  the  materials  conducting  the  propellant. 
We  assume  the  propellant  starts  at  T;0  =  50  K,  and  the  walls,  turbine  blades,  etc.,  limit  the  propellant  high 
temperature  to  Thi  =  3000  K.  This  corresponds  to  an  enthalpy  difference  of  (hhi  —  hi„)  =  105  MJ/kg.  Thus 
we  find  the  mass  flow  of  propellant  by 


Pth 

mprop  hhi  -  hl0 
69.8  MJ/s 
—  105  MJ/kg 

=  0.665  kg/s 


(127) 

(128) 
(129) 


We  divert  the  heated  propellant  to  a  turbine-generator  system,  which  charges  the  capacitor  banks.  If  we 
assume  a  combined  efficiency  of  turbines  and  generators  of  T?jen»7«ur6  =  0.5,  the  power  of  the  propellant  after 
the  turbines  is 


(130) 

Tlgenfhurh 

,  10  MW 

(131) 

69.8  MW - — — 

0.5 

49.8  MW 

(132) 
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We  know  the  propellant  mass  flow  rate  and  the  energy  flow  rate.  We  can  do  a  quick  check  to  see  if  the 
propellant  is  ionized.  The  energy  per  atom  is  roughly  the  propellant  power  divided  by  the  mass  flow  rate: 


kTr 


prop 


as 


prop 

(133) 

TTlprop 

(49.8  MW)  (1.67  x  10"27  kg/atom) 

(134) 

(0.665  kg/s)  (1.602  x  10~19  J/eV) 

0.781  eV 

(135) 

Since  hydrogen  has  almost  no  ionization  at  n  >  1019  m-3  and  kT  as  1  eV,  we  can  use  a  conventional  nozzle 
instead  of  a  magnetic  nozzle.  Thus  we  do  not  need  to  cool  the  ohmic  heating  of  a  magnetic  nozzle’s  coils, 
and  thus  we  can  neglect  Pno*  ■ 

Using  nozzle  characteristics  that  compare  propellant  properties  at  the  turbine  exit,  nozzle  throat,  and 
nozzle  exit, 


kTpr0p 

Xtherm  —  JgT  ^ 

(136) 

_  Vex 

Xv  el  — 

Vthr 

(137) 

we  can  And  the  thrust  and  specific  impulse  of  the  rocket.  Taking 

Xtherm  —  1.35 

(138) 

Xvel  =  2 

(139) 

we  get  kTthr  —  0.578  eV.  The  throat  converts  the  energy  increase  to  kinetic  energy: 

kTthr  kTprop  — 

771  o 
~2Vihr 

(140) 

vthr  = 

2(0.781  -  0.578)  eV(1.602  x  10"19  J/eV) 

1.67  x  10-27  kg/atom 

(141) 

Vthr  = 

6.12  km/s 

(142) 

The  exit  velocity  is 

v„  =  (2)(7.65  km/s) 

(143) 

=  12.2  km/s 

(144) 

The  thrust  and  specific  impulse  are  then 

F 

=  Tfl  propvex 

(145) 

=  (0.665  kg/s)(12.2  km/s) 

(146) 

=  8.14  kN 

(147) 

Up 

Vex 

g 

(148) 

12.2  km/s 

9.8m/s2 

(149) 

=  1250  s 

(150) 

Note  that  while  Up  is  independent  of  v,  F  is  directly  proportional  to  u. 

We  estimate  the  mass  of  the  ship  to  be  the  sum  of  the  payload  and  the  capacitor  banks: 


m, 

mcap 


mpayloa d  "b  mcap 

w 

Peap 


(151) 

(152) 
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where  we  take  pcap  =  2  kJ /kg  to  be  the  specific  mass  of  the  capacitor  banks.  Taking  Tnpaytoad  =  105kg,  we 
get 


105  kg  + 


10  MJ 

2  kJ/kg 


1.05  x  10s  kg 


(153) 

(154) 


If  we  assume  a  tank  fraction  of  Xtank  =0.1,  the  maximum  possible  velocity  change  is  At;ma*  =  29.4 
km/s.  If  we  choose  a  velocity  change,  we  get  the  firing  time  from  Eq.  (24): 

Av  =  25  km/s  (155) 

tjire  =  37.3  days  (156) 

From  Eqs.  (21)  and  (20),  we  can  get  the  initial  rocket  mass,  mo,  and  jet  specific  power,  a j.  We  can  also 
calculate  the  payload  mass  fraction  at  launch: 


a 


2.46  x  106  kg 
4.06% 

49.4  kg/kW 


(157) 

(158) 

(159) 


Pinch  Inventory 


Figure  13:  DPF  Q  vs.  Pinch  Inventory,  N. 

Figure  13  compares  DPF  Q  =  Ef/W  values  for  different  bank  energies,  W,  as  a  function  of  the  number 
of  fusion  fuel  atoms  in  the  pinch,  N.  For  W  =  1  MJ,  the  the  DPF  fusion  relative  output  peaks  around 
N  =  1019,  but  below  breakeven  (Q  <  1).  At  high  pinch  inventories,  the  current  compression  gives  low 
initial  pinch  plasma  temperatures,  ifeTo,  so  the  fusion  reaction  rate  parameter,  <  <rv  >,  is  low.  At  low  pinch 
inventories,  the  current  compression  gives  ifeTo  beyond  the  fusion  <  av  >  peak,  so  the  fusion  rate  is  again 
lower.  At  the  inventory  corresponding  to  the  Q  peak,  the  average  plasma  temperature  during  the  pinch 
lifetime  is  at  the  <  av  >  peak. 

The  peak  Q  increases  as  the  inventory  and  the  bank  energy  increase,  but  for  two  reasons: 

1.  The  increased  bank  energy  gives  a  higher  maximum  current,  Im,  which  overheats  the  initial  plasma 
temperature,  ifeTo-  If  the  inventory  is  increased,  the  average  plasma  temperature  can  be  brought  down 
to  the  <  av  >  peak. 
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2.  Equations  (74)  and  (84)  show  that  the  pinch  lifetime  scales  as  W 1^3.  If  the  pinch  plasma  burns  longer 
at  the  peak  fusion  rate,  the  Q  increases. 

For  W  =  lOA/J.the  DPF  achieves  Q  =  5.30  at  N  =  6  x  1019.  At  higher  bank  energies,  the  peak  Q  levels 
off:  the  pinch  lifetime  is  long  enough,  and  the  fusion  reaction  rate  high  enough,  that  all  the  fusion  fuel  burns. 
Using  the  following  input  criteria, 


w 

— 

10  MJ 

(160) 

io 

= 

1  MA//ts 

(161) 

L 

= 

20  mO 

(162) 

N 

= 

6  x  1019 

(163) 

re 

8  cm 

(164) 

ra 

= 

3  cm 

(165) 

Xtank 

0.1 

(166) 

Figures  14,  15,  and  16  compare  rocket  performance  as  a  function  of  the  firing  frequency,  v,  of  the  DPF.4  The 
hollow  symbols  indicate  the  lighter  payload  mass,  rrifayioad  =  5  x  103  kg  (roughly  the  mass  of  a  satellite), 
and  the  filled  symbols  indicate  the  heavier  payload,  mpavioad  =  105  kg  (heavy  payload).  The  circles  indicate 
a  small  firing  velocity  change,  Av  =  10  km/s  (roughly  the  Av  to  get  to  the  moon),  and  the  triangles  indicate 
a  high  velocity  change,  Av  =  25  km/s  (roughly  the  Av  to  get  to  Mars). 


<U 

6 

£ 
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c 
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Figure  14:  Total  Firing  time  vs.  Firing  Frequency,  v 


As  one  would  expect,  the  low  mass,  low  velocity  change  missions  have  the  shortest  firing  times,  lowest 
jet  specific  powers,  and  the  highest  thrust-to-weight  ratios. 


4  For  systems  where  more  than  one  DPF  fires  at  the  same  time,  the  firing  frequency  would  be  the  product  of  the  number  of 
DPFs  and  the  firing  frequency  of  one  DPF. 


34 


(kg/kW) 


Mpay(kg) 

5e3  le5 
O  •  10  t%r 
V  T  25  (km/s) 


10°  101  102  TO3 


Firing  Frequency  (Hz) 


Figure  15:  Initial  Thrust  to  Weight  Ratio  vs.  Firing  Frequency,  v 
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Figure  16:  Jet  Specific  Power,  a,  vs.  Firing  Frequency,  v 
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A  FORTRAN  source  code 


A.l  heps. for 


C23466789  0 

PROGRAM  haps 

IICLUDE  comb lk. for 


C  begin  main  program 

WRITB(a.*) 'High  Energy  Propulsion  System  (HEPS)  0.99a  06/92* 
URlTE(e.e) 'R.lachtrieb:  OL-AC  PL/RKFE;UIUC  Fusion  Studies  Lab.* 
WRITE(*,e) 'Reactor  models:  DPF.FRC’ 

WRITE(*,*) 

C  read  in  parameters  from  .hep  data  files 
CALL  getvar 

WRITE (*,*) 'Incrementing  variable:  '.IncVar 

WRITE(*, 193) '  lower, upper, howmany  ’ .lower, .upper, ',' .howmany 

*****  begii  calculatiois  ***** 

C  preloop  calculations  conserve  computing  time  by  reducing  work 
C  to  do  in  main  loop. 

CALL  preloop 

IF(loglin.EQ.l)THEl 
C  increment  IncVar  L0G10 

IF(howmany.EQ.l)  THEI 

inc=L0G10(upper)-L0G10 (lower) 

ELSE 

inc= (LOG 10 (upper ) -LOG 10 ( lower) ) / (howmany- 1 ) 

EIDIF 

IF(IncVar . EQ. 'T* )operate*L0G10 (lower) 

IF( IncVar . EQ . ’ deltav ’ ) oper at e=L0G 10 ( lower ) 

ELSE  ! otherwise  increment  IncVar  linearly 

IF  (howmany. EQ.l)  THEI 
inc=upper-lower 
ELSE 

inc* (upper- lower ) / (howaany-1 ) 

EIDIF 

IF(IncVar.EQ. ’T’)operate=lower 
IF ( IncVar . EQ . ' deltav ' ) operat  e=lower 
EIDIF 

**********  BEGII  LOOP  ********** 

DO  40  count* 1, howmany ,1 
bomb*0 

IFdoglin.EQ.  1)THEI  (increment  IncVar  L0G10 

IF(IncVar.Et). ’T’)T*10**operate 
IF ( IncVar . EQ . ’ deltav ' ) delt av* 1 0**oper at e 
ELSE  (otherwise  increment  IncVar  linearly 
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IF ( Inc Var . EQ . ' T  * ) T=oper at  e 
IF ( IncVar . EQ . ' deltav ’ ) deltav=operat  e 
EIDIF 

Ti=T/(i+rt)  ! [keV] 

Te=rt*Ti 


C  plasma  poser  balance,  thermal  posers  lor  reactors 

IF (r eactype . EQ . 2) THE! 

CALL  pos.d 
CALL  therm.d 

c  ELSEIF(reactype.EQ.3)THEI 

c  CALL  pos_l 

c  CALL  therm.f 

EIDIF 

C  system  poser  recirculation,  cooling, 

CALL  cooling 

C  rocket  performance 

CALL  rocket 
IF (bomb . IE . 0) GOTO  39 

C  reactor,  system  masses 

CALL  mass 

*****  EID  CALCULATIONS  ***** 

C  display  output 

CALL  output 

C  goto  here  if  have  to  skip  iteration  (robust) 

39  COITIIUE 

oper at e=oper at  e+ inc 

40  COITIIUE 

**********  EID  LOOP  ********** 

EID 

C  main  program  roc 

**************************** 


BLOCK  DATA 

c  define  constants  in  accordance  sith  F0RTRAI77  syntax 
C  constants 

REAL  muO , k , g , pi , J_MeV , J.keV , u , elec , Uion .NVelec , remnf lux , 

ArhoLiH , NVLiH , Vddn , Vddp , Vd3He , Vpl IB , Wp6Li , Wdt , W23He , Vtt 

C  constants  (from  DATA  blocks) 

CONNOI/COIST/muO , k , g , pi , J.NeV , J.keV , u , elec , Uion , NVelec , remnf lux , 
ArhoLiH , NVLiH , Vddn , Vddp , Vd3He , Vpl IB , VpBLi , Vdt , V23He , Vtt 
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o  o 


C  constants 

DATA  muO , k , g , pi/ 1 . 2566E-6 , 1 . 38066E-23 , 9 . 8062 , 3 . 141692654/ 

DATA  J.MeV . J_keV,u/l . 60219E-13 , 1 . 60219E-16 . 1 . 66056E-27/ 

DATA  elec ,Uion ,NVelsc/l . 602E-19, 13 .6,9. 109534E-31/  ! [C,eV, kg/elec] 
DATA  r emnf lux , rhoLiH , NVLiH/3 . Se-8 , . 82e-3 , 1 . 332e-26/ 

C  [rem  cm '2]  1 or  2.5Mev  neutrons,  like  from  D(d,n)3He 

[kg/ cm* 3] 

MVLiH= [kg/molecule]  molecular  mass  of  lithium  hydride  shield 
DATA  Vddn,Hddp,Wd3He,¥pllB/3.27,4.03, 18.3,8. 68/  ! [MeV/rxn] 

DATA  Vp6Li , Wdt , U23He , Wtt  /4. 02, IT. 69, 12. 9, 11. 3/ 

EID 

************************************************************************* 
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A. 2  comments. for 


C  CONNECTS :  NAZI  VARIABLE  DESCRIPTIONS  (DEPENDENT  AND  INDEPENDENT) 

C  program  incrementable  variables 
C  count=generic  counter 

C  T=Ti+Te 

C  Nr=ratio  propellant  mass  flow  to  fuel  mass  Hob 
C  inc,  increment  f or  DO  loop 

C  expon  [unitless] ,  number  to  which  base  is  raised 
C  functions 

C  SVd3He=[cm*-3  s*-l]  <ev>  parameter  for  D-3He 
C  SVddn= [cm‘-3  s*-l]  <ev>  parameter  for  D(d,n)3He 
C  SVddp= [cm*-3  s*-l]  <ev>  parameter  for  D(d,p)T 

C  constants 

C  muO=1.2666E-6;  [H/m]  permeability  of  free  space 

C  k=1.38068E-23;  CJ/K]  Boltzman’s  constant 

C  J_keV=l . 60219E-16 ;  [J/keV] 

C  u=l .66056E-27;  [kg]  atomic  mass  unit 

C  g  *  9.8062;  [m/s‘2]  gravitational  constant 

C  pi-3.141592654;  c 

C  elec=l .602E-19  [C]  electron  charge 

C  HWelec=[kg]  electron  mass 

C  rhoLiH=[kg/cm“3]  density  of  lithium  hydride  (shield  material) 

C  HWLiH=[kg]  molecular  mass  of  lithium  hydride  shield 

C  remnf lux- [rem  cm'2]  for  2.5N eV  neutrons,  like  from  D(d,n)3He  (w2.45NeV) 
C  Wd3He=[J]  energy  produced  per  D(3He,p)'  fusion  reaction 
C  Vddn=[J]  energy  produced  per  D(d,n)3He  fusion  reaction 
C  Wddp=[J]  energy  produced  per  D(d,p)T  fusion  reaction 
C  Uion=[eVj .ionization  potential  (13.6  for  hydrogen) 


C  INDEPENDENT  VARIABLES 
C  reactor 

C  reactype=  type  of  reactor  (DPF.FRC,  etc) 

C  lower  [keV]  lower  temperature  boundry  for  DO  loop 

C  upper  [keV]  upper  temperature  boundry  for  DO  loop 

C  howmany=number  of  iterations  to  perform 
C  Zj  [unitless] .atomic  number  of  element  j 
C  Zk  [unitless] .atomic  number  of  element  k 
C  Aj  [unitless] ,  atomic  mass  of  element  j 

C  Ak  [unitless] ,  atomic  mass  of  element  k 

C  kj  [unitless] .fractional  density  of  element  j 
C  kk  [unitless] .fractional  density  of  element  k 
C  rt  [unitless] ,  ratio  T(electron)/T(ion) 

C  Kc  [unitless] .fraction  of  cyclotron  radiation  absorbed  in  Beil Is 
C  cycabs  [unitless] .fraction  of  cyclotron  radiation  absorbed  in  plasma 


46 


C  y=fraction  of  charged  particle  power  retained  in  plasma 
C  Xn=ratio  o f  propellant  density  to  fuel  density  in  reaction  chamber 
C  etaGen=generator  efficiency  [Wth  ->  We] 

C  rl®[cm]  inner  shield  radius  (4c  type  shield) 

C  rpssngr=[cm]  passenger  distance  from  neutron  source 
C  Xcone=fraction  of  sphere  conical  shield  occupies 
C  specRad [kg/Wth] =  specific  mass  of  radiators 
C  specGen [kg/Wel] =  specific  mass  of  generators 

C  specExt Dcg/Wel] =  specific  mass  of  external  electrical  eng.  supply 
C  resmat=[jm]  resitivity  of  theta-pinch  material 
C  rhomat3 [kg/cm‘3]  density  of  theta-pinch  material 
C  tensmax3 [Pa]  tensile  strength  of  theta-pinch  material 

C  Xtens=ratio  of  theta  pinch  material  tensile  strength/stress  (safety  factor) 
C  maxdose=[rem]  maximum  dose  allowable  for  person  from  reactor 
C  ddnspin-suppressing/ enhancing  factor  from  spinpolarization  (l=no  effect) 

C  ddpspin3suppr easing/ enhancing  factor  from  spinpolarization 
C  d3Hespin-suppressing/ enhancing  factor  from  spinpolarization 
C  etaInj=injector  efficiency  [We  ->  Wplasma] 

C  speclnj Ckg/Wmdf uel] =  specific  mass  of  fusion  fuel  injectors 

C  dpf 

C  PICHRAD=[cm]  radius  of  pinch 
C  RAIODE=AIODE  RADIUS  (cN) 

C  RCATH=CATHODE  RADIUS  (cN) 

C  LAIODE3 ABODE  LE1GTH  (cM) 

C  PICHRAD=RADIUS  OF  PIICH  (cM) 

C  LPICH  =PIHCH  LEIGTH  (cM) 

C  FS*PLW3SIOWPLOW  efficieicy  factor,  fractioi  of  iiitial  fill  gas 
C  WHICH  IS  EITRAIHED  II  THE  RUIDOWI. 

C  FPICH  =PIICH  EFFICIEICY  FACTOR.  FRACTIOI  OF  GAS  II  RUIDOWI  WHICH 
C  IS  TRAPPED  IISIDE  THE  PIICH. 

C  RHOIOI=IIITIAL  FILL  GAS  DEISITY  (KG/cm**3) 

C  VOLT3CHARGIIG  POTEITIAL  OF  CAPACITOR  BAIKS  (V) 

C  CAP= IIITIAL  EXTER1AL  CAPACITAXCE  (CAPACITOR  BARK)  (F) 

C  LIlITsIIITIAL  IIDUCTAICE  OF  EXTERIAL  CIRCUIT  (H) 

C  ITERS3IUMBER  of  iteratiois  performed  DURIIG  EACH  PIICH 
C  REPRATE-IUMBER  OF  FIRIIGS  PER  UIIT  TIME  (S**-l) 

C  PICHTIM=DURATIOI  OF  STABLE  PIICH  PHASE  (S) 

C  DSCBRG=TIME  FOR  FILL  GAS  TO  BE  DISCHARGED  (S) 

C  IMAGIET=CURREHT  IECESSARY  TO  PRODUCE  2  TESLA  MAGIETIC  FIELD  (A) 

C  specCap=SPECIFIC  EIERGY  OF  CAPACITOR  BAIKS  (J/KG) 

C  volmix=[cm*3]  volume  of  mixing  chamber  (for  DPF  only) 

C  volmix=volcham  for  FRC 
C  Bnoz3CT]  nozzle  magnetic  field 

C  frc 

C  choice=particle  confinement  scaling  law  choice  (Krall,  VSLS,  TRX,  etc) 

C  dzplug*[cm]  thickness  of  mirror  opposing  nozzle-end  mirror. 

C  mult*fraction  of  gl(VSLS)  achievable 
C  rate=parameter  for  finding  gl  for  VSLS  theory 
C  near spar ameter  for  finding  gl  for  VSLS  theory 
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C  Lcham= [cm]  chamber  length 
C  Is  Cca] ,  separatrix  length 

C  xs  [unitless] .ratio  separatrix  radius/coil  radius 
C  rcoil  [cm] ,  radius  of  FRC  coil 

C  XBplug=ratio  plug  magtetic  field  strength  to  theta  pinch  fid  strength 
C  XBnoz=ratio  nozzle  magtetic  field  strength  to  theta  pinch  fid  strength 
C  Xrplug=ratio  inner  radius  of  plug  to  coil  radius 
C  Xmoz=ratio  inner  radius  of  nozzle  to  coil  radius 

C  rocket 

C  Zprop= [unitless] , atomic  number  of  propellant 
C  Aprop=*  [unitless] ,  atomic  mass  of  propellant 
C  hO=[J/kg]  initial  enthalpy  of  cold  propellant 

C  hMat=[J/kg]  enthalpy  of  propellant  at  high  temperature  (materials  limited) 
C  Xtemp=Tmix/Tthroat  p  Xtemp  w  1.3S  i.e.  is  nozzle  characteristic 
C  Xvel=VEXIT/VTHROAT  p  Xvel  s  2.0  i.e.  is  nozzle  characteristic 
C  tmd=[days]  mission  time 
C  payload= [kg]  payload  mass 

C  shelter3 [kg]  mass  of  Solar  Particle  Event  (SPE)  and  storm  shelter 
C  Ppayld*[W]  poser  demands  of  payload 
C  dznoz=[cm]  thickness  of  nozzle  coil 

C  output 
C  oO  > 

C  ol  > 

C  o2  >  I/O  toggles  to  specify  shat  to  srite  to  file 
C  .  > 

C  .  > 

C  olO  > 


C  for  dpf.  n  is  dependent 

C  Bnoz  is  dependent  for  FKC 

C  volmix  is  dep.  for  FRC 


C  DEPEXDEIT  VARIABLES  (passed  to  subroutines) 

C  plasma  parameters 
C  Ti  [keV] ,  average  ion  temperature 
C  Te  QceV] .  average  electron  temperature 
C  volplas=[cm"3]  plasma  volume 

C  NVfuelsQcg]  molecular  mass  of  average  fuel  particle 
C  fz  [unitless],  =kj*Zj+kk*Zk ,  describes  number  of  electrons  in  plasma 
C  n  [cm* -3] .total  ion  density  in  plasma 
C  mdf uel* [kg/ s]  fuel  mass  flos  rate 

C  beta  [unitless] ,  kinetic  plasma  press/external  magnetic  field  press 
C  Bint=[T]  internal  magnetic  field  strength 
C  Bext=[T]  external  magnetic  field  strength 
C  RRd3He=[cm*-3  s*-l]  reaction  rate  of  D(3He,p)' 

C  RRddn=[cm*-3  s*-l]  reaction  rate  of  D(d,n)3He 
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C  RRddp*[cm*-3  s*-l]  reaction  rate  of  D(d,p)T 
C  taul*[s]  confinement  time  of  particles  in  FRC  plasma 

C  frc  parameters  (dependent) 

C  kappa*parameter  for  finding  gl  for  VSLS  theory 
C  drcoil=[cm]  theta  pinch  coil  thickness 
C  rplug= [cm]  plug  inner  radius 
C  drplug3 [cm]  plug  thickness 
C  Bpluga [T]  plug  magnetic  field  strength 

C  dpf  parameters  (dependent) 

C  IMAXSQ*MAXIMUM  CURREIT  SQUARED  (AMPS) 

C  IHAX=sqrt  MAXIMUM  CURREIT  SQUARED  (AMPS) 

C  I MOPT* ASSUMED  VALUE  OF  MAXIMUM  ATTAIIABLE  CURREIT  (AMPS) 

C  TOTEFF=TOTAL  EFFICIEHCY=FSIPLW*FPHCH 
C  XSECTAR*CROSS  SECTIONAL  AREA  OF  FOCUS  DEVICE  (cM**2) 

C  IFTHRST=THRUST  DUE  TO  EXPELLED  (ION-PINCH)  GASES  (I) 

C  ELECTEI =T0T AL  ELECTRICAL  ENERGY  TO  BE  SUPPLIED  BY  CAPACITORS  (V) 

C  powers 

C  Preactor* [W]  necessary  power  to  drive  reactor 

C  Pinj*[W]  power  injected  into  plasma  (if  nec.  to  maintain  steady  state) 

C  Pbr  [W] ,  bremstrahhlung  radiation 
C  Pcy  [V] ,  cyclotron  radiation 

C  Pradwall=[W]  thermal  radiation  power  absorbed  (to  be  dissipated) 

C  Pth=[W]  total  thermal  power  to  be  dissipated 
C  PelReq= [W]  electrical  power  required 
C  PelGen* [W]  power  produced  by  generators 

C  PradRf 1= [W]  cyclotron  rad  reflected  (abs.  by  prop  and  plasma) 

C  Pradplas=[W]  power  absorbed  by  plasma  from  PradRf 1 
C  Ppayload=[tf]  power  requirements  of  payload 

C  Pcoil*[W]  power  necessary  to  produce  Bcoil  (also  nec.  to  dissipate  I*2yR) 
C  Pprop* [W]  power  of  heated  propellant  flow. 

C  Pmix* [W]  thermal  power  at  mixnation  point 
C  Pip3 [V]  power  contained  in  particle  lost  from  plasma 
C  PextReq*[W]  external  electrical  power  required 

C  Pnoz=[W]  power  necessary  to  produce  Bnoz  (also  nec.  to  dissipate  I~2yR) 

C  Pplug=[V]  power  necessary  to  produce  Bplug  (also  nec.  to  dissipate  I*2yR) 
C  TFP=TOTAL  FUSION  POWER  (W) 

C  u=ratio  thermal  radiation  absorbed  by  propellant/total  thermal  rad. 

C  gamGen*f r act ion  of  heated  propellant  sent  to  generators 
C  gamCool-fraction  of  cold  propellant  sent  to  absorb  thermal  rad. 

C  masses 

C  reactor3 [kg]  reactor  mass 
C  injector* [kg]  mass  of  fuel  injector 
C  powersys=[kg]  mass  of  power  system 
C  MCAP=TOTAL  CAPACITOR  BANK  MASS  (KG) 

C  prop* [kg]  total  propellant  mass 
C  fuel* [kg]  total  fuel  mass 
C  nozzle* [kg]  nozzle  mass 
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C  radiator- Deg]  radiator  mass 
C  shi«ld=  Deg]  shield  nass 
C  mf = [kg]  final  system  mass 
C  mO=Deg]  initial  system  sums  =  mf+prop+fuel 

C  rochet  parameters  (dependent) 

C  raoz=[cm]  magnetic  nozzle  inner  radius 
C  drnoz=[cm]  magnetic  nozzle  thickness 
C  Tmix*Dt]  nozzle  miznation  temperature 

C  n_mix» [cm*-3]  density  of  particles  at  miznation  (for  nozzle) 

C  thrusts [I]  total  thrust  produced 
C  Isp=[s]  specific  impulse 
C  tms=[sec]  mission  time 

C  HWprops  Deg]  molecular  mass  of  propellant  particle 
C  mdprop= Ckg/s]  propellant  mass  flow  rate 
C  mdout=  Dcg/s]  total  mass  flow  rate  out 
C  ionized=fraction  of  ionized  plasma 
C  Bnoz=[T]  nozzle  magnetic  strength 

C  component 

C  surfsratio  surface  area  avail,  cyclo.  rad.  /total  chamber  surface  area 
C  volchams Ccm'3]  reaction  chamber  volume 
C  rl= [cm]  distance  of  shield  from  reactor 
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A.3  comblk.for 


IMPLICIT  REAL  (a  -  z) 

CHARACTER  choice*l , IncVar*8 , out_name*12 ,pause*l 

LOGICAL  00,01,02,03,04,06,06,07,08,09,010,011,012,013,014,015 

DOUBLE  PRECISIOI  ionized 


C  connon  blocks  sharsd  by  subroutines 

C  progran  increnentable  variables 

COMMOI/PROG/ count , IncVar , hownany , lover , upper , inc , operate , bomb , 
kloglin 

C  constants  (Iron  DATA  blocks) 

COMMOI/COIST/muO , k , g , pi , J.MeV , J.keV , u , elec , U ion , MVelec , rennf lux , 
ArhoLiH .MWLiH , Wddn , Vddp , Vd3He , Wpl IB , Vp6Li , Udt , U23He , Vtt 
C  independent  variables  (read  in  iron  .hep  Tiles) 

C QMM01/IKD EP/oO, o 1, o2, o3, o4, 06, 06, o7, o8,o9,ol0,o ll,o 12,o 13, ol4, 
to 15 , choice , mult , rat  e , near , Lchaa , Is , r co il , zs , ZBplug , XBnoz , Xrplug , 
kdzplug , RAIODE , RCATH , LA10DE , RHOIOI , VOLT , CAP , LIIIT , FSIPLV , FPICH , 
ALPICH , PICHRAD , REPRATE , PICHTIM , ITERS , DSCHRG , INAGIET , specCap , 
Rvolnix , r eactype , rt , Kc , cy cabs , y , In , etaGen , r 1 , rps sngr , Xcone , 
AspecRad , specGen , specExt , resnat , rhonat , t  ensnax , Xt  ens , naxdose , 
Addnspin,ddpspin,d3Hespin,apecInj , etalnj ,Zprop,Aprop,hO,hMat, 
AXtenp , Xvel , Ppayld , dznoz , Xraoz , payload , shelter , Xtank , T , deltav , 
Rxprot , xdeut , xtrit , x3He , x6Li , xl IB , rail a 
C  dependent  variables  (passed  to  subroutines) 

CQHMQl/DEPEI/Ti ,Te , volplas , KVfuel , Iz ,n , ndf uel , beta, Bext .Bint , 
ARRddn , kappa , dr coil , rplug , drplug , Bplug , taul , 

ATOTEFF , XSECTAR , IFTHRST , IMAX , ELECTEI , P  j  etIF , 

APreact or , Pin j , Phr , Pcy , PradRf 1 , Pradplas , Pradwall , Pth , Pip, Pnix , 
APelGen , PelReq , PextReq, Pnoz , Pcoil , Pplug , TFP , v , ganCool , ganGen , 

Ain j  ector , Neap , reactor , radiator .shield , nozzle , poser sy s ,  nf , nO , 
Aprop ,  fuel ,  moz ,  dmoz ,  Bnoz ,  Tnix ,  n.nix ,  thrust ,  Isp ,  ionized, 

Atna , MVprop .ndprop , ndout , surf , volchan , tanks 


C  for  iterations  tracking 

91  FORMAT  (A.F6.0.A) 

92  FORMAT  (A,F5.0,A,A,A,1PE11.3) 


C 

C 

131 

132 

133 

134 

135 

136 

137 

138 


10  TEXT  FORMAT 

three  places  after  decimal 
FORMAT  (1PE11.3) 

FORMAT  ( 1PE1 1.3, 1PE11 . 3) 

FORMAT  (1PE11 .3,1PE11.3,1PE11.3) 

FORMAT  ( 1PE1 1 . 3 , 1PE1 1 . 3 , 1PE1 1.3.1PE11.3) 

FORMAT  (1PE11.3,1PE11.3,1PE11.3,1PE11.3,1PE11.3) 

FORMAT  (1PE11.3,1PE11.3,1PE11.3,1PE11.3,1PE11.3,1PE11.3) 

FORMAT  (1PE11.3,1PE11.3,1PB11.3,1PE11.3,1PE11.3,1PE11.3,1PE11.3) 
FORMAT  (1PE11.3,1PE11.3,1PE11.3,1PE11.3,1PE11.3,1PE11.3,1PE11.3, 
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A1PE11.3) 

139  FORMAT  (lPE11.3,lPE11.3,lPE11.3,lPE11.3.1PE11.3.1PE11.3tlPE11.3, 
A1PE11 . 3, 1PE11 . 3) 

130  FORMAT  (1PE11.3,1PE11.3,1PE11.3,1PE11.3,1PE11.3,1PE11.3.1PE11.3, 

A1PE11 . 3, 1PE11 . 3 , 1PE11.3) 

C  TEXT  FORMAT 

191  FORMAT  (A.1PE11.3) 

192  FORMAT  (A,1PE11.3,A,1PE11.3) 

193  FORMAT  (A,1PE11.3,A,1PE11.3,A,1PE11.3) 
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A.4  getvar.for 


C234S6789  0 

SUBROUTIIE  getvar 

IKCLUDE  conblk.for 

*****  IIPUTS  ***** 

C  read  in  data  Iron  files 

C  45  reactor 

□PEI  (UIIT*46,  FILE* 'REACTOR. hep* ,  STATUS* ’ OLD  * ) 

READ (45,*) 

READ (46 , * ) reactype 
READ (46,*) 

READ(45,'(A8)')IncVar 
READ (46,*) 

READ  (  45 ,  *  )  hosnany 
READ(46,*) 

READ (46,*) lover 
READ(45,*) 

READ (45,*) upper 
READ(46,«) 

READ (46 , * ) logl in 
RSAD(45,*) 

READ(46,*)T 

READ(46,*) 

READ(46,*)xprot 

READ(46,*) 

READ(46,*)xdeut 
READ (46,*) 

READ(45,*)ztrit 
READ (46,*) 

READ(45,*)z3He 

READ(46,*) 

READ(46,*)x6Li 
READ(46, *) 

READ(45,*)zllB 
READ(45, *) 

READ(45,*)zalfa 
READ  (45,*) 

READ(46,*)rt 

READ(46,«) 

READ(46,*)Kc 

READ(45,*) 

READ (45 , * ) cycabs 
READ (45,*) 

READ(45,*)y 

READ(45,*) 

READ(46,*)Xn 
READ  (46,*) 

RBAD(45,*)*taGen 
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READ  (45,*) 

READ(45,*)rl 
READ (46,*) 

READ ( 45 , * ) rps s&gr 
READ (45,*) 

READ (45 , * ) Xcona 
READ  (45,*) 

READ(45 , *) spacRad 
READ (46,*) 

READ (45 , * ) apacGan 
READ (45,*) 

READ (45 , * ) apacExt 
RBAD(45,*) 

READ (46 , * ) raaaat 
READ  (45,*) 

READ (45 , * ) rhomat 
READ (45,*) 

READ (46 , *) tanamax 
READ (46, *) 

READ(45,*)Ztana 
READ (45,*) 

READ(45 , *)aaxdoaa 
READ (46.*) 

READ (46 , * ) ddnapin 
READ  (45,*) 

READ(46, *)ddpapin 
READ(46,*) 

READ (46, *)d3Haapin 
CLOSE  (UIIT=46) 

C  41  frc 

IF(raactypa . EQ . 3) THEY 

□PEI  (UIIT=41,  FILE3 'FRC . hap ’ ,  STATUS3 'OLD') 
READ (41,*) 

RBAD(41, ’ (Al) ' )choica 
RBAD(41,*) 

READ(41,*)nralt 

READ(41,*) 

READ(41,*)rata 
READ (41,*) 

READ(41,*)naar 

READ(41,«) 

RBAD(41,*)Lcha» 

READ(41,«) 

READ(41,*)la 

READ(41,*) 

READ(41,*)rcoil 
READ (41,*) 

READ(41,*)xa 

READ(41,*) 

READ(41,*)n 
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READ (41,*) 

READ (41 , *)XBplug 
READ (41,*) 

READ(41,*)XBnoz 
READ (41,*) 

READ(41,*)Xrplug 

READ(41,*) 

READ(41,*)dzplug 
READ (41,*) 

READ(41 , *)Xrnoz 
READ (41,*) 

READ (41 , * ) spec In j 
READ(41,*) 

READ(41,*)etaInj 
CLOSE  (UIIT=41) 

E1DIF 

C  43  dpi 

IF(r*actype . EQ . 2) THE* 

OPE*  (U*IT=43 , FILE* ' DPF . Hap  * , STATUS* 'OLD') 
READ (43,*) 

READ(43,*)RAI0DE 

READ(43,*) 

READ ( 43 , * ) RC ATH 
READ  (43,*) 

C  kr 

C  no  lpnch 

READ (43 , * ) LAIODE 
READ (43,*) 

READ(43,*)RH0I0I 
READ (43,*) 

READ (43,*) VOLT 
READ(43,*) 

READ  (43.*)  CAP 
RBAD(43,*) 

READ(43,*)LIIIT 

READ(43,*) 

READ (43 , * ) FSIPLH 
READ (43,*) 

RBAD(43,*)FPICH 

READ(43,*) 

READ(43,*)LPICH 

READ(43,*) 

READ(43,*)P*CHRAD 

READ(43,«) 

READ (43,*) REPRATE 
READ(43,*) 

READ(43,*)PICHTIN 

READ(43,*) 

READ (43,*) ITERS 
READ(43,*) 
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READ (43 , * ) DSCHRG 
READ (43,*) 
READ(43,*)IMAGIET 
READ(43,*) 
READ(43, *)Bnoz 
READ (43,*) 

READ (43 , * ) specCap 
READ (43,*) 
READ(43,*)Xmoz 
READ (43,*) 

READ (43,*) apecln j 
READ(43,*) 

READ (43 , * ) etaln j 
CLOSE  (U1IT=43) 
BID  IF 


C  47  rocket 

OPE*  (U1IT=47 ,  FILE=' ROCKET. hep' ,  STATUS= ' OLD ’ ) 
READ(47,*) 

READ (47, *) Mr 
READ  (47,*) 

RBAD(47,*)Zprop 

READ(47,*) 

READ (47,*) Apr op 
READ(47,*) 

READ(47,*)hO 
READ  (47,*) 

READ(47,*)hMat 

READ(47,*) 

READ(47,*)Xtenp 

READ(47,*) 

READ(47,*)Ival 

READ(47,*) 

READ(47,*)deltav 
READ (47,*) 

READ ( 47 ,*) payload 
RBAD(47,*) 

RBAD(47,*)Ppayld 

READ(47,*) 

READ(47,*) shelter 
READ(47,*) 

READ(47,*)dznoz 

READ(47,*) 

READ (47 , * ) volaix 
READ(47,*) 

READ(47,*)Ztaak 

READ(47,*) 

READ(47,*)Tkix 
CLOSE  (UIIT*47) 


C  49  output 
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OPEI  (U*IT=49,  FILE=» OUTPUT. hep 

READ (49, 

*) 

READ (49, 

*)o0 

READ (49, 

*) 

R£AD(49, 

*)ol 

READ (ew. 

*) 

READ (49, 

*)o2 

READ (49, 

*) 

READ  (49, 

*)o3 

READ  (49, 

*) 

RBAD(49, 

*)o4 

READ (49, 

*) 

READ(49, 

*)o6 

READ (49, 

*) 

READ(49, 

*)o6 

READ (49, 

*) 

READ (49, 

*)o7 

READ (49, 

*) 

READ (49, 

*)o8 

READ (49, 

*) 

READ (49, 

*)o9 

READ  (49, 

*) 

READ(49, 

*)olO 

READ (49, 

*) 

READ (49. 

*)oll 

READ (49, 

*) 

READ (49, 

*)ol2 

READ (49, 

*) 

READ (49, 

*)ol3 

READ(49, 

*) 

READ (49, 

*)ol4 

READ(49, 

*) 

READ (49, 

*)ol6 

CLOSE  (UIIT=49) 

STATUS= ’ OLD  * ) 


*****  DISPLAY  IIPUTS  ***** 

C  display  input  values  to  check  file  read  was  accurate 
IF(oO)THBI 

WRITE(*,*) 'Press  <Eater>  to  continue’ 

READ(*, ’ (Al) * ) pause 

C  46  reactor 

WRITE(* , *) ’Data  fro*  file:  REACTOR’ 

VRITEC*,*) ’Incrementing  variable:  ’.IncVar 

VRITB(*,191) ’  howmany  * ’ , howmany 

WRITE(*,191) ’  lower  lower 

WRITE (* ,191)’  upper  * * , upper 

WRITE(*, 191) ’  T  CkeV]  «’,T 

WRITE(*,191) ’  zprot  *’,xprot 

WRITE(*, 191) ’  /leut  *’,xdeut 
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WRITER,  191)’  xtrit  =’,xtrit 

WRITEC*, 191)’  x3He  =’,x3He 

WRITEC*, 191) ’  x6Li  =’,x6Li 

WRITEC*, 191) ’  xllB  =  ’ ,xllB 

WRITEC*, 191) ’  xalla  =’,xalia 

WRITEC*, 191) ’  rt  =’ ,rt 

WRITEC*, 191) ’  Kc  =' ,Kc 

WRTTEC*, 191) ’  cycaba  cycaba 

WRITEC*, 191) *  y  ,y 

WRITEC*, 191) ’  In  =  ’,Xn 

WRITEC* , 191) '  etaGen  ■’.etaGen 

WRITEC*, 191)'  rlCca]  =',rl 

WRITEC*, 191) '  rpaangrCca]  =’,rpsangr 

WRITEC*, 191)'  Icon*  =’,Xcone 

WRITEC*, 191) '  specRad Qcg/kWth]  =’,apecRad 

WRITEC*, 191) ’  apecGen [kg/kWel]  =',specGen 

WRITEC*, 191)'  apecExt [kg/kWel]  =', apecExt 

WRITEC*,*) 'Preaa  <Enter>  to  continue' 

REAOC*, ' (Ai) ' )pauae 

WRITEC*,*) 'Data  Iron  file:  REACTOR  (continued)' 
WRITEC*, 191) '  reamatCjm]*’ .reamat 
WRITEC*, 191) ’  rhomat [kg/ cm*  3] * ' .rhomat 
WRITEC*, 191) '  tenamax[Pa]=’ .tenamax 
WRITEC*, 191)'  Xtena  =',Xtena 

WRITEC*,  191)’  maxdoaeCrea]  =’,maxdoae 
WRITEC*, 191) ’  ddnapin  =' .ddnspin 
WRITEC*, 191)’  ddpapin  *’ .ddpapin 
WRITEC*, 191)'  d3Heapin  *',d3Hesp in 
WRITEC*,*) 'Preaa  <Ent*r>  to  continue’ 

REAOC*, ’ (Al) ’ ) pause 

C  41  Ire 

IFCreactype . EQ • 3)THE> 

WRITEC*,*) 'Data  *rom  lile:  FRC’ 

WRITEC*,*) 'choice  choice 

WRITEC*, 191)’  mult  »’,mult 

WRITEC*, 191)’  rate  =’,rate 

WRITEC*, 191) ’  near  =' .near 

WRITEC*, 191)’  LchamCcm]  * ’ , Lcham 

WRITEC*, 191)’  la  Ccm]  *',la 

WRITEC*, 191)’  r coil [cm]  =’,rcoil 

WRITEC*, 191)'  xa  =  ’,xa 

WRITEC*, 191)’  n[cm*-3]  *’,n 

WRITEC  * ,191)’  XBplug  * ' . XBplug 

WRITEC* , 191) ’  XBnoz  * ’ . XBnoz 

WRITEC*, 191)’  Xrplug  * ’ , Xrplug 

WRITEC*, 191)'  dzplug Ccm]  = ' , dzplug 

WRITEC*, 191) ’  Xrnoz  * ’ , Iraoz 

WRITEC*, 191) ’  apeclnj [kg/kWmdfuel]  *’ .speclnj 

WRITEC*, 191)’  etalnj  *’,etalnj 

WRITEC*,*) 'Preaa  <Enter>  to  continue’ 
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READC*.  'CA1)  ’)pause 

EIDIF 

C  43  dpf 

IFCreactype . EQ . 2)THEM 
WRITEC*,*) 'Data  Iron  iila:  DPF' 

WRITEC* , 191) '  RAMODE  [ca] =’, RAMODE 
WRITE ( *,191)'  RCATH  Cca] = ' , RCATH 
WRITE(*, 191) '  LAIODE  [cm]=' .LAIQDE 
WRITE(*, 191) '  RHOIOV  =’,RHOIOM 
WRITEC*. 191) ’  VOLT  [V] = ' , VOLT 
WRITE(*,191) ’  CAP  [F]»',CAP 
WRITEC*. 191) ’  LIMIT  [H]  =’ .LIMIT 
WRITEC*, 191) ’  FSMPLW  =’,FSIPLW 
WRITE(* , 191) '  FPMCH  =',FPSCH 
WRITEC*, 191) ’  LPMCH  [cm]=',LPMCH 
WRITE(* , 191) ’  PMCHRAD  [ca]=’ .PHCHRAD 
WRITEC*. 191)'  REPRATE  Cs*-1]=' .REPRATE 
WRITEC*, 191)'  PMCHTIM  Ca]=’ .PICHTIM 
WRITEC*, 191) ’  ITERS  ='. ITERS 
WRITEC*, 191)'  DSCHRG  [a]=’,DSCHRG 
WRITEC*. 191)'  IMAGHET  =« .IMAGHET 
WRITEC* ,*) 'Press  <Enter>  to  continue' 

READC*,  ’  CA1)  ') pause 

WRITEC*,*) 'Data  from  file:  DPF  C continued) ’ 

WRITEC*, 191)’  Bnoz  [T]=’,Bnoz 

WRITEC*, 191) ’  specCap  [K j /kg] = ' , apecCap 

WRITEC*.  191)  ’  Xrnoz  =’, Xrnoz 

WRITEC*. 191) ’  speclnj  =’,speclnj 

WRITEC*, 191)'  etalnj  =',etalnj 

WRITEC*,*) 'Press  <Enter>  to  continue' 

READC*, ’ CA1) ' )pause 

EMDIF 

C  47  rocket 

WRITEC*,*) 'Data  from  file:  ROCKET’ 

WRITEC*. 191)’  Hr  =’,Hr 

WRITEC *.191)’  Zprop  =’,Zprop 

WRITEC*, 191) ’  Aprop  =’,Aprop 

WRITEC*. 191)’  hO  [J/kg]  =  ’,h0 

WRITEC*, 191)’  hMatCJ/kg]  *’ ,hMat 

WRITEC*. 191) ’  Xteap  =', Xteap 

WRITEC*. 191)’  Xvel  =’,X7el 

WRITEC*, 191) ’  deltavCa/a]  =’,deltav 

WRITEC*. 191) ’  payloadCkg]  =', payload 

WRITEC*. 191)’  PpayldCW]  *',Ppayld 

WRITEC*, 191) ’  shelter [kg]  =’, shelter 

WRITEC*. 191)’  dznoz  [cm]=’,dznoz 

WRITEC*, 191) ’  volaix  Ccm*3]=' .volmix 

WRITEC*. 191)’  Xtank  =’,Xtank 

WRITEC*, 191) ’  TmixtK]  Cinitial  gueas)=’ .Trnix 
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WRITE(*,*) ’Press  <Enter>  to  continue' 
READ(*,'(Al)’)pause 

EIDIF 

E1D 

C  end  subroutine  getvar 
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A. 5  preloop.for 


C23456789  0 

SUBROUTIHE  preloop 

IICLUDE  comb lk . for 

*****  BEGII  GEIERIC  PRE-LOOP(EXT)  CALCULATIOIS  ***** 

mui=xprot*l . +xdeut *2 . +xtrit*3 . +x3He*3 . +x6Li*6 . +xl IB* 1 1 . +xalf a*4 . 

C  [vmitleea]  average  atomic  mass 

MWfuel=u*mui 

C  [kg]  ave  mass  of  tael  ion 

f z=xprot+xdeut+rtrit+x3He*2 . +x6Li*3 . +xl 1B*5 . +xalf a*2 . 

C  fz=(kj*Zj+kk*Zk+. . .) 

C  [unitless]  average  atomic  no  o f  particles 

KWprop=u* (Apr op) 

C  [kg]  mass  o 1  propellant  ion 

*****  EID  GEIERIC  PRE-LOOP (EXT)  CALCULATIOIS  ***** 

WRITE(*,*)  ’lame  of  file  to  output?  (8+3  char,  nu)  ' 

READ(*, ’ (A12) ’)out_name 

□PEI  (UIIT=65,  FILE=out_name ,  STATUS='*EW') 

WRITE(*,*) 'Writing  to  file  '.out.name 


IF(reactype . EQ . 2)THEI 
WRITE(56,*)’DPF’ 

*****  BEGII  DPF  PRE-LOOP  CALCULATIOIS  ***** 

C  CALCULATE  FOCUS  PARAMETERS 
lpnch  =  rcath  -  r anode 
pnchrad  =  ranode/kr 

XSECTAR=pi*(RCATH**2 .  -  RAI0DE**2.) 

C  [cm“2] 

volcham=lanode*XSECTAR 
C  [cm*3] 

volplas=pi*PICHRAD**2 . *LPICH 
C  [cm"  3] 

TOTEFF =FSIPLW  *  FPICH 

surf  = ( 2 . *pi*RCATH*LAIODE  +  2 . *pi*RAHODE*LAHODE  +  pi*RCATH**2 . ) 
k/  (2 . *pi*RCATH*LAIODE  +  2 . *pi*RAIODE*LAIODE  +  2.*pi*RCATH**2. ) 
C  fraction  of  surface  area  not  "holes",  ie  ends 

*****  EID  DPF  PRE-LOOP  CALCULATIOIS  ***** 

ELSEIF (reactype . EQ . 3)THEI 
WRITERS,*)  ’FRC’ 

*****  BEGII  FRC  PRE-LOOP (EXT)  CALCULATIOIS  ***** 
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c 


beta=l-xs**2./2. 

FRC  formula 

volplas=ls*pi* (xs*rcoil) **2 . 

C  cm  "3 

volcham=Lcham*pi*rcoil**2 . 

C  cm*  3 

C  compute  particle  confinement  time  using  one  of  three  scaling  methods. 
WRITE(*,*) 'Particle  Confinement  Scaling:  ’ 

IF  (choice. Eq. 'V' .OR. choice. Eq. 'v')  THE1 
C  VSLS  Theory 

WRITE(S5,*)'VSLS  Theory' 

WRITE(*,*)’VSLS  Theory' 
kappa=BETA 

10  IF  ((TAIH(kappa) /kappa). GT. BETA)  THEI 

kappa=kappa*RATE 
ELSE 

kappa=kappa/RATE 

EIDIF 

DIFF=ABS (BETA  -  TAIH (kappa) /kappa) /beta 
IF  (DIFF . GT . BEAR)  GOTO  10 
ELSEIF  (choice. Eq.'T'. OR. choice. EQ.'t')  THE! 

WRITE(56,*)'TRX-2  Experiment' 

WRITE(*,*) 'TRX-2  Experiment' 

ELSEIF  (choice. Eq. 'K*. OR. choice. Eq. 'k')  THEI 
WRITE(55,*) 'Krall  Theory* 

WRITE(*,*)'Krall  Theory' 

ELSEIF  (choice. Eq.'L'. OR. choice. Eq.'l*)  THEI 
URITE(S5 , * ) ' Loser  Hybrid  Drift  Theory' 

WRITE(*,*) 'Lower  Hybrid  Drift  Theory' 

EIDIF 

*****  EID  FRC  PRE-LOOP  (EXT)  CALCULATION  ***** 

EIDIF 

C  write  header  to  table-format  output  file 

IF(o2.AID.o9)VRITE(68,*)'T,  qppPf us/Pinj ' 

IF(o3. AID. o9)WRITE(65.*)'T,mdout, prop, fuel' 

IF(o4. AID.o9)URITE(56,*) 'T, injector, radiator, shield, reactor, ’ , 

A ' nozzle , poser sys , payload , shelt er ' 

IF ( 06 . AID . o9 ) WRITE ( 56 , * ) ' T , Pr adsall , P co il , Pnoz ' 

IF(o6 . AID . o9) WRITE (55 , *) 'T, Ppayld , Pin j /etaln j , Pcoil , Pplug, Pnoz ' 
IF(o7. AID.o9)WRITE(55,*) 'T,s,gamCool,gamGen,PelGen/PextReq' 
IF(o8.AID.o9)WRITE(5E,*) 'T,  Thrust,  Ion-Fusion  Thrust,  ratio' 

c  IF (o 10. AID. o9) WRITE (66,*) 

c  ft ' t ime , xprot , xdeut , xtr it , x3He , x6Li , xl IB , xalf  a , n ’ 

c  IF(olO. AID.o9)WRITE(55,*) 

c.  ft ' t ime , PFddn , PFddp , PFd3He . PFpl IB , PFpflLi , PFdt , PF23He , Pneut , PFTOT ' 

IF(olO . AID . o9)WRITE(65 ,*) 'T,TFP, TPneut ,TFP/TPneut ' 


IF(ol2. AID. o9) WRITE (55,*) 
ft' T,  deltaV,  tmd,  payload/mO' 
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IF(ol3.AID.o9)WRITE(55,*) 

*'T,  mdluel*Vrun**2/2,  TDELTAP ,  Pprop,  Pmix,  Pjet' 

IF(ol4. AID .o9)WRITE(65,*) 

*'T,  P j etlF/Pmix ,  Pc/Pmix,  Pprop/Pmix,  Pjet[W]' 

IF(ol5. AID.o9)VRITE(5S,*) 

*’T,  Pjet,  PjetIF,  Pinj' 

c  IF(olO)THEI 

c  VRITE(*,*)  ’ VARITOP  destination  file  name?  (8+3  char,  max)  ' 

c  READ(*, ' (A12) ' )out_na»e 

c  OPEI  (UIIT=65,  FILE=out .name ,  STATOS= ’ IEW ' ) 

c  WRITE(*,*) 'Writing  to  file  *,out_nane 

c  EVDIF 


EKD 

C  subroutine  preloop 
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A.6  gettau.for 

SUB ROUTI IE  gettau 
C  finds  tauX 

C  truly  local  variables 

REAL  const , Fkrall , R , omegaci , rhoi , Scap , stau , Vtaa , f  s 

IMPLICIT  REAL  (a  -  z) 

CHARACTER  choice*l , IncVax*8 

LOGICAL  00,01,02,03,04,05,06,07,08,09,010,011,012,013,014,015 
DOUBLE  PRECISIO!  ionized 

IICLUDE  comblk.for 


*********  pRC  PARTICLE  COVFIIEMEIT  SCALIIG  taul  *************** 

IF  (choice. Eq. 'T' .OR. choice. Eq. 't')  THEI 
C  TRX-2  Experiment 

const=l.e-10 

C  calibrates  scaling  lam  to  give  confinement  times  in  proper  units 

taul= (xs*rcoil)**l . *(n/l . el5) **0 . 9*T**1 . 5*const 

from  L.Steinhauer  personal  conversations,  Tnsz.  FRC  review 
a/lelB  converts  n  to  [10~15  cm"-3]  for  scaling  lav 
ELSEIF  (choice. Eq. *V» .OR. choice. Eq. *v')  THEI 
C  VSLS  Theory 

const =6. e-6 

C  calibrates  scaling  lav  to  give  confinement  times  in  proper  units 

taul= (Ti* 1 . e3) **0 . 2*kappa**2 . 7  *rcoil** 1 . 8*Bext**2 . 4* 
ft  (n/l.el6)**(-l.)*const*mult 

C  from  Miley  luc.Instm.Meth. 

C  Ti*le3  converts  keV  to  eV 

C  n/lel6  converts  n[cm*-3]  to  n[10*15  cm~-3]  for  scaling  lav 

ELSEIF  (choice. Eq. 'K' .OR. choice. Eq. 'kO  THE! 

C  Krall  Theory 

Fkrall =pi/ (3 . * ( 1 . +rt ) **0 . 5 ) 

taul=3 . e-8*(Bext*l . e4)*(xs*rcoil)**2. /(Fkrall*(Te*l . e3)* 
ft  (l.+rt)**0.6) 

C  from  Krall  Phys. Fluids. 

C  input  taulCs],  BextCT],  xs*rcoil[cm] ,  T[keV] 

C  Ti*le3  converts  keV  to  eV 

C  Bext*le4  converts  T  to  Gauss  for  Krall' 8  scaling 

ELSEIF  (choice. Eq.'L'. OR. choice. Eq.'l’)  THEI 
C  Lover  Hybrid  Drift  Theory 

R*(rcoil*xs/100. )/(2.**0.5) 
omegaci=elec*Bext/MWfuel 
rhoi=(J_keV*Ti/MWfuel)**0.5/omegaci 
Scap=R/rhoi 

stau*xs**2 . *Scap/2 . **(3 . /2 . ) *3 . /4 . 

C  corresponding  to  delO  from  Phys  Fluids  26(1983)1626  Hoffmann/Milroy 

Utau=3 . 
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C  between  2-4,  but  taul  not  very  sensitive  to  Vtau  anyway 

fs=l. 

C  "arbitrary"  parameter 

taul=0 . 13*f s** (-1 . ) ebeta* (stau**3 . ) *Scap/( 1 . +Te/Ti)/Bint* 

*  (l.+Wtau/(stau/3.))**3./l.e6 

C  le6  converts  taufsCfs]  to  tanlCs] 

C  from  luc  Fus  24(1984)1540  Slough  et  al.  eq.  7  for  spatially 

C  varying  (LHD)  resistivity  profile 

E1DIF 

C  now  have  taul 
EID 

C  real  function  taul 


A.7  sigmav.for 


C  given  ion  teaperature  (keV) ,  escalates  sigma-v  parameter 
Coefficients  thanks  to:  Larry  T.  Cox,  AL(AFSC)/LSVF 

C  for  D(d,n)3He  rxn; 

REAL  FUXCTIOI  SVddn(T) 

REAL  a,b,c,d,T,x,L0G10 
x*L0G10(T) 
a*  0.2981119 
b*  -2.082958 
c=  5.701349 
d=-22. 08780 

SVddn= 10** (a*x**3+b*x**2+c*x+d) 

BID 

C  function  SVddn 


C  for  D(d,p)T  rxn; 

REAL  FU1CTI01  SVddp(T) 

REAL  a,b,c,d,e,T,x,L0G10 
x=L0G10(T) 
a=  -0.05730514 
b=  0.5716992 
c=  -2.404634 
d=  5.643689 
e=-21. 97106 

SVddp=10** ( a*x**4+b*x**3+c*x**2+d*x+e) 
BID 

C  function  SVddp 


C  for  3He(d,p)‘  rxn; 

REAL  FUICTIOI  SVd3He(T) 

REAL  a,b,c,d,T.x,L0G10 
x=L0G10(T) 
a=  0.3535589 
b=  -3.310354 
c=  10.10471 
d=-25. 67344 

SVd3Ha='10**(a*x**3+b*x**2+c,*x+d) 

BID 

C  function  SVd3He 

C  for  ilB(p.‘)'+‘  rxn; 

REAL  FUICTIOI  SVpllB(T) 

REAL  a.b.c.d.T.x.LOGlO 
x=L0G10(T) 
a=  0.6879438 
b*  -6.833501 
c*  17.40498 
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d=-33. 51638 

SVpUB=10**(a*x**3+b*x**2+c*x+d) 

EKO 

C  {unction  SVpllB 

C  lor  6Li(p,3He)‘  rxn; 

REAL  FUVCTZOI  SVp6Li(T) 

REAL  a,b,c,d,T,x, LOG 10 
x=L0G10(T) 
a=  0.5712202 
b»  -4.265839 
c»  11.87649 
d=-28. 12494 

SVp6Li=10**(a*x**3+b*x**2+c*x+d) 

E1D 

C  lunction  SVp6Li 

C  lor  T(d,n)'  rxn; 

REAL  FUICTIOI  SVdt(T) 

REAL  a.b,c.d,T.x,L0G10 
x=L0G10(T) 
a^  0.2708006 
b=  -2.421732 
c=  6.318869 
d=-20. 15761 

SVdt=10**(a*x**3+b*x**2+c*x+d) 

E1D 

C  lunction  SVdt 

C  lor  3He(3He,2p)'  rxn; 

REAL  FUICTIOI  SV23He(T) 

REAL  a,b,c.d.T,x,L0G10 
x=L0G10(T) 
a=  0.6392849 
b=  -5.274763 
c=  16.43308 
d=-l 1.10836 

SV23He=10**(a*x**3+b*x**2+c*x+d-24) 

c  -24  above  is  lor  barns 

EID 

C  lunction  SV23He 

C  lor  T(t,2n) ‘  rxn; 

REAL  FUICTIOI  SVtt(T) 

REAL  a,b,c,d,e,T,x,L0G10 
x*L0G10(T) 
a*  -0.07238949 
b=  0.6854508 
c=  -2.827826 
d=  6.513504 
e»-22. 72402 
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SVtt=10**(a*x**4+b*x**3+c*x**2+d*x+e) 

BID 

C  function  SVtt 
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A.8  pow-d.for 


C23466789  0 

SUBROUTIIE  pow.d 

IICLUDE  comblk.for 

*****  BEGII  DPF  PLASMA  POWER  BALAICE  ***** 

C  PLASMA  TEMPERATURE 

AU*(W*IdotO*Ldot*Ldot)**(i./3. ) 

AI= ( W*IdotO/Ldot ) ** ( 1 . /3 . ) 

At= (W/ ( IdotO*IdotO*Ldot ) ) ** ( 1 . /3 . ) 

I max  =  0.64  *  AI 

voltO  =  2.12  *  AU 

voltmax=  0.64  *  AU 

tris*  =  1.46  *  At 

C  RU1D0WI  VELOCITY 

C  THE  RUIDOWI  VELOCITY  IS  CALCULATED  USIIG  THE  MOMEITUM  EQUATIOI 
Vrun=Ldot*2*pi/ (muO*log(rcath/ranod«) )  ! Cm/ a] 

mdlu*l=HWfuel*inven*REPRATE  ! Ckg/a] 

C  DETERMIIE  IIJECTED  POWER  BY  RATE  OF  BAIK  El ERG I ES  DELIVERED 
Pinj=W*REPRATE  ! [W] 

C  initial  loop  values 
C  PIICH  HUMBER  DEISITY 

n=inven/volplas  ! [cm*-3] 

C  define  initial  element  number  densities  in  pinch  via  specified  fractions 
prot=n*xprot 
deut=n*xdeut 
trit=n*xtrit 
hel3*n*x3Be 
lith=*n*x6Li 
boron=n*xllB 
alfa=n*xalfa  ! [cm* -3] 

f z=xprot+xdeut+xtrit+x3He*2 . +x8Li*3 . +xl 1B*6 . +xalf a*2 . 
fzSQ=xprot+xdeut+xtrit+x3He*4. +x6Li*9. +xllB*26 . +xalf a*4. 

TFP=0 

TPL0SS=0 

TDELTAP=0 

Pbr*0 

Pcy*0 

TPn*ut=0 

iprot  =xprot 
ideut  =xdeut 
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itrit  =xtrit 
i3He  =x3He 
i6Li  =x6Li 
illB  sxl IB 
ialla  =xalla 

C  FOB  OIE  PULSE 

*****  begin  DPF  internal  loop  ***** 

DO  2002  count2=l, iters, 1 

IF (MOD ( 100 . *count2/iters , 5 . ) . EQ . 0 . ) 

*  WRITE(*,91) '+  dpi’  ,count2/iter«*100. , 

C  to  determine  progress 

RRddn  =deut*(deut  *SVddn(T)  ) 

RRddp  =deut*(deut  *SVddp(T)  ) 

RRd3He=deat* (hel3  *SVd3He (T) ) 

RBpllB=prot*(boron*3VpllB(T) ) 

RRp6Li=prot*(litb  *SVp6Li(T) ) 

RRdt  =deut*(trit  *SVdt(T)  ) 

RR23He=hel3* (hel3  *SV23He(T)) 

RRtt  =trit*(trit  *SVtt(T)  )  !  [cm*-3  s*-l] 

vrpi=volplas*REPRATE*PICHTIM/ITERS 

C  DETERMIIE  CHARGED  FUSIOI  POWER  FROM  PIICH 
PFddn  =vrpi*J_MeV*RRddn*Wddn*0.2S 
PFddp  =vrpi*J_MeV*RRddp*Wddp 
PFd3He= vrpi*  J_MeV  *RRd3He*Wd3He 
PFpllB=vrpi*J_MeV*RRpllB*WpllB 
PFp6Li=vrpi*J_MeV*RRp6Li*WpflLi 
PFdt  =>vrpi*J_MeV*RRdt*Wdt*0.20 
PF23He=vrpi*J_MeV*RR23He*W23He 
PFtt  =vrpi*J_MeV*RRtt*Wtt 

Pneut  =vrpi* J_MeV* (RRddn*Wddn*0 . 75+RRdt*Wdt*0 . 80) 
PFT0T=PFddn+PFddp+PFd3He+PFpllB+PFp6Li+PFdt+PF23He+PFtt  ! [W] 

C  average  atomic  number  ol  particles  in  pinch 


fz=iprot+ideut+itrit+i3He*2 . +i6Li*3 , +il 1B*5 . +ialf a*2 . 
lzSq=iprot+ideut+itrit+i3He*4 . +i6Li*9 . +illB*2S . +ialla*4 . 

PBREM=n*(n*6 . 3SE-31)*lz*lzSq*SqRT(Te) *vrpi 
C  n[c*‘-3] ;  TeCkeV];  Pbrem.Pcyc  [W] ; 

prot=prot+(-RRpllB-RRp6Li+RRddp+RRd3He+2.*RR23He)*PICHTIM/ITERS 
deut=deut+ (-RRdt-2 . * RRddn- 2 . *RRddp-RRd3He)*PICHTIM/ITERS 
trit=trit+(-RRdt-2 . *RRtt+RRddp) *PICHTIM/ ITERS 
h*13=hel3+ ( -RRd3He-2 . *RR23He+RRddn+RRp6Li) *PICHTIM/ITERS 
lith=lith+ ( -RRpflLi) *PICHTIM/ ITERS 
bor on=boron+ ( -RRp 1 IB ) *PICHTIM/ ITERS 
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alia*alfa+(+RRdt+RRtt+RRd3He+RRp6Li+3.*RRpllB+RR23He) 
t  *PICHTIM/ITERS 

n=prot+deut+trit+h«13+,litli+boroa+aIla 

iprot  aprot/n 

ideut  =deut/n 

itrit  =trit/n 

i3He  =hal3/n 

i6Li  =lith/n 

illB  =boron/n 

ialfa  aalfa/n 

TFP=TFP+PFTOT 
TPneut=TPneut+Pneut 
Pbr=Pbr+PBREM 
Pcy=Pcy+PCYC 
C  CW] 

2002  COVTIIUE 

*****  and  internal  loop  ***** 

IF ( o 10 ) WRITE ( 55 , 134) T , TFP , TPneut , TFP/TPneut 

TDELTAP=TFP-Pbr 
Plp=TDELTAP 
C  CW] 

*****  EHD  DPF  PLASMA  POWER  BALAICE  ***** 

BID 

C  subroutine  pow_d 
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A.9  therm-d.for 


C23466789  0 

SUBROUTIHE  th«rm_d 

IICLUDE  comblk.for 

*****  BEGII  DPF  THERMAL  POWER  ***** 

C  CYCLOTRO*  RADIATIQI  ABSORBED  II  WALL  AID  ELECTRODES 
Pradsall=Pbr+Kc*Pcy 

C  umat  sails  ars  r*fl*ctiv*  to  Bremsstrahlung  radiation 

C  TOTAL  POWER  DISSIPATED  II  ELECTRODES  BY  OHMIC  HEATIIG 
Pelctrd=Imax**2. *resmat*(LAI0DE/100 . )*(1 . / (pi*( 

ARAIODE/ 100 . ) **2 . ) 

*+l . / (pi* ( (RCATH/100 . ) **2 . -(RAI0DE/100 . ) **2 . ) ) ) *DSCHRG*REPRATE 
C  TOTAL  POWER  TO  BE  REMOVED  FROM  THE  WALLS  AID  ELECTRODES  AID  MAGIET 

c  on  first  run,  Tmix  is  not  yet  calculated:  naad  nozzla  thermal 

c  thermal  poser  to  determine  Tmix,  and  Tmix  to  determine 

c  nozzle  thermal  poser 

c  check  Tmix,  ionization;  if  out  of  bounds,  say  so;  if  need 

c  conventional  nozzle,  Bnoz  =  0,  drnoz  fixed  at  XX 

c  drnoz  minimum  =  1  cm 

C  nozzle  strength  Bnoz  depends  on  mixnation  temperature  and  ion 

C  density  of  effluent 

C  Bnoz- 10 . **(0 . 5*L0G10(Tmix [eV] )+L0G10(n_mix)-17 . 05) 

C  from  "Char.  Flos  Thru  Mag.  Ioz“ ,  Sgro,  Glasser,  et  al. 

C  n_ions3n_fuel  +  *n_prop 

n_ ions= (n*volplas+Xn*n* ( volcham-volplas ) ) / volcham 
C  [cm* -3] 

C  n_elec*n_fuel*Z_fuel  +  n_prop*Z_prop 

n_elec= ( n*volplas*f z+Xn*n* ( volcham-volplas ) *Zprop) / volcham 
C  [cm* -3] 

n_mix= (n_ions+n_elec) *volplas/ volmix 

Bnoz= 10 . ♦* (0 . 5*L0G 10 (Tmix/ 1 1604 . 86 ) +L0G 10 (n_mix) -17 . 05) 

C  calculate  DPF  magnetic  nozzle  strength  for  file  input. 

C  don't  calculate  above  because  read  in  Bnoz  (Tmix,n_mix 

C  depend  on  Pnoz,  Bnoz  i.e.  not  calculated  yet 

C  ultimately  only  affects  mass 

C  USE  SAHA  EqUATIOI  TO  CALCULATE  DEGREE  OF  PR0PELLAIT  I0IIZATI0I 
C  FOR  MAGIETIC  I0ZZLE  (assumes  hydrogen  propellant) 

C  if  ionization  isn't  near  100%  mag.  nozzle  is  not  as  effective 
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IF  (Tmix.LT.2S00)  THEN 
C  use  conventional  nozzle 

draoz=10 . 

C  [cm] 

Pnoz=0 . 

C  CW] 

ELSE 

CALL  s ah a  (Uion,Tmix,n_mix, ionized) 

IF  ( ionized. LT. 0.3)  THE! 
c  WRITE(* , *) *  tf  AS_D :  *Hot  ionized  ' , 

c  *’ enough  for  magnetic  nozzle*' 

drnoz=20 . 

C  [cm] 

Pnoz=0 . 

C  CW 

ELSE 

rnoz=Xmoz*rcath 

draoz=moz/(tensmax/Xtens*2.  *muO/Bnoz**2.-l ./3. ) 

IF  (dmoz.LT.l)  drnoz=l 
alphg=  (rnoz+dmoz)  /moz 
betag=dznoz/ (2*moz ) 

gab= (betag/ (2 . *pi* (alphg**2 . -1 . ) ) ) **0 . 5* 

4L0G ( ( alphg+ ( alphg**2 . +betag**2 . ) **0 . 5 )/( 1 . + ( 1 . +betag**2 . ) **0 .5)) 
lambda=0 . 9 

Pnoz=(Bnoz/ (lambda*muO*gab) ) **2 . *reamat*moz/100 . 

C  rnoz/100.  converts  moz  [cm]  to  raozCm] 

ENDIF 

ENDIF 

Preactor=Pelctrd 

Pth=Pradwall+Preactor+Pnoz 

*****  END  DPF  THERMAL  POWER  ***** 

END 

C  subroutine  was_d 
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A. 10  pow-f.for 


C23468789  0 

SUBROUTIHE  pou_f 

COMHOI/. . . 

*****  BEC.II  FRC  PUSH!  POWER  BALAICE  ***** 

C  plasma  pressure  balance:  Press_mag,ext=Press_plas  +  Press.mag, int 

C  and  def  <a>pPress_plas/Press_mag,ext=l.-xs**2./2. 

Bext=(n*l . e6* J_keV*T/(beta/(2 . *muO) ) ) **0 . S 
C  n*l.e6  converts  n[cm*-3]  to  n[m*-3] 

Bint=((Bext**2./(2.*muO)-n*l .e6*J_keV*T)*2.*mu0)**0.5 

C  find  particle  confinement  time 

C  tauI(l8,xs,rcoil,n, const, rs=xs*rcoil,B,Te,Ti, kappa, etc. . .) 

CALL  gettau 

C  D(3He ,p) '  D(d,n)3He  D(d,p)T 

RRddn=ddnspin*(Kj *n) **2 . /2 . *SVddn(T) evolplas 
RRddp=ddpspin* (K j *n) **2 . /2 . *SVddp(T)*volplas 
RRd3He=d3Hespin*kj*kk*n**2 . *SVd3He(T)*volplas 
C  [cm" -3  s‘-l] 

Pc=RRd3He*Wd3He+RRddp*Wddp+RRddn*Wddn/4 . 
TFP=RRd3He*Wd3He+RRddp*Wddp+RRddn*Wddn 
C  [W] 

C  RRddn*Wddn/4 .  is  fraction  of  DDn  power  to  3He,  not  n 

C  steady  state,  fuel  losses  equal  inputs 

C  RRd3Be . . .  can  be  considered  negligible  compared  to  n/taul,  but  include 

C  to  be  accurate 

escrate  =n*volplas/tauI 
C  [s*-l] 

bumrate=RRd3He+RRddn+RRddp 
C  [s"-l] 

lo8srate=escrate+bumrate 
C  [s*-l] 

mdfuel=lo8srate*MWfuel 
C  [kg/s] 

surf=(2.*pi*rcoil*Lcham  +  pi*(rcoil**2.-(Xrplug*rcoil)**2. ) 
k+  pi*(rcoil**2.-(Xraoz*rcoil)**2. )) 

A/(2. *pi*rcoil*Lcham  +  2. *pi*rcoil**2. ) 

C  fraction  of  surface  area  not  "holes"  (ends) 

C  calculations  of  radiation  sources 

cl=6. 35e-31*fz* (k j  *Z j **2 . +kk*Zk**2 . ) 

dl=2 . 5e-32*f z**2 . * (1 . +Ti/ (f z*Te) )  *<1 . +Te/204 . ) 

Pbr=cl*n**2 . *Te**0 . 6*volplas 
C  [W] 
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Pcy=dl* ( 1 . -beta) /beta*n**2 . *Te**2 . *Kc*volplas 
C  CM] 

C  assume  that  all  of  reflected  cyclotron  radiation  is  absorbed  by  fuel 
C  plasma  and  propellant :  assume  it  is  divided  between  fuel  and  prop 
C  according  to  H 
ref 1=1. -Kc 

PradRfl=refl*surf*Pcy 
C  [W] 

Pradplas=cycabs*PradRfl 
C  [V] 

C  plasma  power  balance:  inputs=outputs  (steady  state) 

C  using  enter/ezit  control  volume 

C  Pinj+Pc=Plp+Pbr+(Pcy-Pradplas) 

Pin j  =Plp+Pbr+ (Pcy-Pradplas ) -Pc 
IF  (Pinj.LT.O)  Pinj=0. 

C  can’t  have  negative  injected  energy  -  excess  energy  will  go  to 

C  heat  up  plasma:  what  happens  to  plama  temperature,  particle  loss  rate 

C  ?Examine.  Can  hold  at  certain  temp  by  varying  fuel  mixture 

Plp=(l.-y)*Pc  +  escrate*J_KeV*T 
C  [W] 

*****  EHD  FRC  PLASMA  POUER  BALA ICE  ***** 

EHD 

C  subroutine  pow_f 
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A. 11  therm-f.for 


SUBROUTIHE  therm_f 
C0MM0I/ . . . 

*****  BEGII  FRC  THERMAL  POWER  ***** 

C  examine  later  if  can  ignore  neutron  heating  of  shield 

C  else  have  to  get  into  KERMAs  (look  at  Steve  Gruneisen  paper) 

Prad=Pbr+Kc*Pcy 

C  CALCULATE  I*2R  HEATXIG  1H  THETA  PIICH,  PLUG,  AID  IOZZLE 
C  Iron  Dolan  p.611,  tensile  stress  on  solenoid® 

C  eT=(B"2/210) (rl/(r2-rl)+l/3) ,  drcoil=r2-rl ,rl=rcoil 
drcoil=rcoil/ (t ensmax/Xtens*2 . *muO/Bext**2 . -i . /3 . ) 
alphg® (rcoil+drcoil)/rcoil 
betag=Lcham/(2. *rcoil) 
gab=(betag/ (2 . *pi*(alphg**2 .-1 . ) ) )**0 . 5* 

ALOG ( ( alphg+ ( alphg**2 . +betag**2 . ) **0 . 5 ) / ( 1 . + ( 1 . +betag**2 . ) **0 . 6) ) 
lambda=0 . 9 

Pcoil= (Bert/ (lambda*muO*gab) ) **2. *resnat*r coil/100. 

C  rcoil/100.  converts  rcoilCcm]  to  rcoil[m] 

C  assume  plug  area  is  square,  thickness  of  plug 
Bplug=XBplug*Bext 
rplug=Xrplug*rcoil 

drplug=rplug/ (t ensmax/Xtens*2 . *muO/Bplug**2 . -1 . /3 . ) 

alphg= (rplug+drplug) /rplug 

betag=dzplug/ (2*rplug) 

gab= (betag/ (2 . *pi* (alphg**2 . -1 . ) ) ) **0 . 5* 

ALOG ( ( alphg* ( alphg**2 . +betag**2 . ) **0 . 5 ) / ( 1 . + ( 1 . +betag*  *2 . ) **0 . 5 ) ) 
lambda-0 . 9 

Pplug® (Bplug/ (lambda*muO*gab) ) **2 . *resnat*rplug/100 . 

C  rplug/ 100.  converts  rplug [cm]  to  rplug Cm] 

Bnoz=XBnoz*Bplug 

rnoz=Xrnoz*rcoil 

dmoz=moz/  ( t  ensmax/Xt  eus*2 .  *muO/Bnoz**2 .  - 1 .  /3 . ) 
alphg=  (rnoz+dmoz)  /rnoz 
betag®dznoz/ (2 . *rnoz) 
srite(*f *)  ’rnoz.draoz’  ,moz,dmoz 
writ  e ( * , * ) ' alphg , betag ' , alphg , betag 
write(* , *) ' (betag/ (2 . *pi*(alphg**2 . -1 . ) ) ) **0 . S ' 
write(*,*)(betag/(2.*pi*(alphg**2.-l.)))**0.5 
gab=(betag/(2.*pi*(alphg**2.-l . )))**0.5* 

ALOG ( (alphg* (alphg* *2 . +betag**2 . ) **0 . 5) / (1 . + ( 1 . +betag**2 . ) **0 . 5) ) 
lambda=0 . 9 

Pnoz= ( Bnoz/ ( 1 ambda*mu0*gab ) ) **2 . *re  smat  *rnoz/ 100 . 

C  rnoz/ 100.  converts  rnoz Com]  to  rnoz  fc] 

Preactor=Pcoil+Pplug 
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Pth=Prad+Preactor+Pnoz 


C 

***** 


C 


[W] 

EHD  FRC  THERMAL  POWER 


EID 

subroutine  sas_f 


***** 
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A. 12  cooling.for 


SUBROUTIIE  cooling 
IICLUDE  comblk . lor 

*****  BEGII  GE8ERZC  COOLIVG,  POWER  GEI. ,  Mix.  COVOS.  ***** 

C  SET  PROPELLAIT  FLOW  VALVES  TO  COOL  COKPOIEITS,  RU1  GEIERATOR 

C  II  sending  all  the  propellant  to  be  heated  absorbs  all  the  thermal 

C  radiation  poser,  only  as  much  propellant  sill  be  sent  that  sill  alios  the 
C  highest  temperature  difference  (materials  limited)  to  make  the  generators 
C  more  efficient. 

C  Similarly,  if  sending  all  heated  propellant  to  generators  produces  too 

C  much  elec,  poser  then  only  as  much  prop,  as  nec.  sill  be  sent; 

C  otherwise  the  remainder  is  dumped  in  cold. 


c«*****«****g0  RADIATOR 

mdprop=Pth/  (hJfat  -hO) 
mdout=mdfuel+mdprop 
C  [kg/s] 

gamCool-1 
s=l 

hmin=0. l*hmat 
gamGen=l 

sdot =gamGen*mdprop* (hMat-hmin) 
PelGen=etaGen*sdot 

PelReq=Ppayld+Preactor+Pnoz+Pinj/etaInj 
C  [W] 

PextReq=PelReq-PelGen 
IF  (PelGen . GT . PelReq)  THE! 
PelGen=PelReq 
sdot-PelGen/etaGen 
ga*G*n=sdot/(mdprop* (hMat-hmin) ) 
PextReq*0 
BID  IF 


C  solve  for  enthalpy  of  propellant  mixing  together  again  after  cooling 
C  nozzle,  running  generator,  or  being  dumped  in  cold 
hprop=hKat * ( 1 . -gamGen)  +  (hMat-hmin) *gamGen 
Pprop=mdprop* (hprop-hO ) 

C  [W] 

C  Assume  the  plasma  loss  particles  (sith  poser  Pip)  and  and  the  propellant 

C  (sith  poser  Pprop)  mix  uniformly  and  come  to  a  mixing  temperature  Tmix 

C  (i.e.  velocity*0). 

Pmix*Plp+Pprop 
C  [W] 
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*****  EUD  GEIERIC  COOLIIG,  POWER  GEM.,  Mix.  CORDS. 


EXD 

C  subroutine  cooling 


***** 


79 


A. 13  saha.for 


SUBROUTIVE  saha  (0, temp, density, ionized) 

DOUBLE  PRECISIOI  dpi ,dp2,dp3,dp4,Xsaha, ionized 
REAL  EXP 

REAL  U.temp, density 


C  need  to  used  double  precision  because  of  limited  word-length  of 

C  IBM  PC  (le+/-37) 

C  n*lD6  converts  n[cm*-3]  to  a [a" -3] 

C  temp/1 1604. 85  converts  TCK]  to  TCeV] 

C  U  is  required  ionization  energy,  13.6  eV  for  hydrogen 

dpl=(U/(temp/11604.85)) 
dp2=(3 . 313e-28) *(density*l . e6) 

IF(temp.LT.O.)  THE* 

write(*,*) *Saha:  Hot  ignited* 
ionized=0 
RETURI 
ELSE 

dp3= (t emp/ 1 1 604 . 85 ) ** ( -3 . /2 . ) 

EIDIF 

IF(dpl .GT.70 . )  THEI 
ionized=0 . 

ELSE 

dp4=EXP(dpl) 

Xsaha=dp2*dp3*dp4 
IF  (Xsaha.EQ.O)  THE! 

vrite(*,*) 'Saha  algorithm  failed...' 
write(*,*) ,dpl3(U/(temp/11604.85))=* ,dpl 
write(*,*) *dp2=(3.313e-28)*(density*l . e6)=* ,dp2 
write(*,*) ’dp3=(temp/11604.85)1»*(-3./2. )=’ ,dp3 
write(*,*) *dp4=EXP(dpl)=* ,dp4 
ionized=-l . 

ELSE 

ionized* ( ( 1 . +4 . *Xsaha) **0 . 5-1 . ) /(2 . *Xsaha) 

EHDIF 

EIDIF 


E1D 

C  subroutine  saha 
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o  o  o 


A. 14  rocket. for 


SUBROUTIIE  rocket 
IICLUDE  comblk . for 

*****  BEGII  GEIERIC  ROCKET  CALCS.  ***** 

C  THRUST  AID  SPECIFIC  IMPULSE 

Assume  propellant  particlea  enter  into  the  throat  of  the  magnetic 
nozzle  completely  ionized;  therefore  the  magnetic  nozzle  can  direct 
the  propellant  along  with  the  plaama  lose  particlea. 

C  n_iona=n_fuel  +  *n_prop 

n_iona= (n*volplaa+Xn*n* (volcham-volplas ) ) / volcham 
C  [cm*-3] 

C  n_elec=n_fuel*Z_fuel  +  n_prop*Z_prop 

n_elec=(n*volplas*fz+Xn*n*(volcham-volplas)*Zprop)/volcham 
C  [cm* -3] 

n_mix=(n_ions+n_elec)*volcham/volmix 

y_ions=n_ions/(n_ions+n_elec) 

y_elec=n_elec/(n_iona+n_elec) 

C  mole  fraction  of  iona , electrona 

C  MWiona=(I_fuel*MVfuel  +  l_prop*MVprop)/I_iona 

MVions= (n*volplae*MWfuel+Xn*n*MVprop* (volcham-volplas) ) 

A  /(n_ions*volcham) 

C  [kg/ ion  (ave)] 

x_iona=y..ione*HWiona/(y_iona*KWion8  +  y_elec*HVelec) 
x_elecay_elec*M¥elec/(y_iona*MWiona  +  y_elec*HVelec) 

C  mass  fraction  of  iona .electrons 

C  non  the  power  is  divided  among  an  electron  gas  and  an  ion  gas. 

C  R=8. 3143  J/mol  K 

C  assuming  an  ideal  gas,  Cp=Cv  *  R 

C  ie.  Cp=R*(5./2.0)=20.786  J/mol  K 

C  but  R/Iav=k=Boltzmann 1  a  Constant 

CpIon=5 . /2 . *k/MWions 
CpElec=5 . /2 . *k/MWelec 

C  ions  and  electrons  come  to  thermal  equilibrium  and  all  energy  is 
C  thermal  (stagnation) 

C  Pmix=(Cp*Tmix)mdout,  solve  for  Tmix 

Tmix=Pmix/(mdout*(CpIon*x_iona  +CpElec*x_elec) ) 

C  Tmix/TTH  p  Xt«p  n  1.36  AID  VEXIT/VTHROAT  p  Ivel  n  2.0 

TTH=Tmix/It«mp 


C  USB  SAHA  EQUATIOI  TO  CALCULATE  DEGREE  OF  PROPELLAYT  I0IIZATI0I 
C  FOR  MAGIETIC  I0ZZLE  (aaaumea  hydrogen  propellant) 
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C  if  ionization  isn’t  n«ar  100%  mag.  nozzle  is  not  as  effective 


CALL  saha  (Uion,Tmix,n_mix, ionized) 

C  from  conservation  of  energy:  CP*DELTAT= ( 1/2) *VTHR0AT**2 
C  thermal  energy  is  converted  to  enthalpy  of  the  plasma 
VionTH=SQRT(2. *CpIon*(Tmix-TTH) ) 

V  el ecTH= SQRT ( 2 . *CpEle  c  * (Tmix-TTH ) ) 

C  floe  exits  Zvel  times  as  fast  because  of  expansion 
C  through  the  nozzle. 

V ionEX=Xvel*V ionTH 
VelecEX=Xvel*VelecTH 

elecTHR=x_elec*mdout*VelecEX 
ionTHR=x_ions*mdout*VionEX 
IF (reactype . EQ . 3 ) IFTHRST=0 . 

C  if  FRC,  no  "non-pinch”  thrust 

thrust=elecTHR  +  ionTHR  +  IFTHRST 
Isp=thrust/ (g*mdout ) 

IF ( Isp*g*L0G ( 1 . /Xtank+1) . LE . deltav) THEI 

VRITE(*,*) 'ROCKET:  Isp  too  lorn  for  requested  delta  V; 
A'Try  higher  energy.’ 
bomb=l 

C  prevent  bomb  in  subsequent  subroutines  with  SQRT(-..) 

EKDIF 

*****  ESD  GEIERIC  ROCKET  CALCS.  ***** 


E1D 

C  subroutine  rocket 
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A. 15  mass.for 


SUBROUTIHE  mass 
IICLUDE  comblk. for 

*****  BEGII  SPECIFIC  REACTOR  MASS  CALC  ***** 

IF(r eactype . EQ . 2)  THEM 
*****  BEGII  DPF  MASS  CALC  ***** 

C  DETERMIVE  THE  IECESSARY  CAPACITOR  MASS  USI1G  IHPUTTED 
C  SPECIFIC  EIERGY 

C  ASSUME  EFFICIE1CY  OF  COIVERSIOI  FROM  ELECTRICAL  TO 
C  MAGIETIC  EIERGY  IS  COISTAIT. 

C  VO-O . 5*C0*V0**2=B**2/ (2*mu0) 

C  SIICE  B  IS  PROPORTIOIAL  TO  I,  ASSUME  IECESSARY  STORED 
C  ELECTRICAL  EIERGY  IS  PROPORTIOIAL  TO  1**2. 

C  H=V0*(I0PT/I0)**2 

C  IEED  ELECTEI  II  kJ  AID  specCap  II  kJ/kg 
ELECTEI=0 . S*CAP*V0LT**2 . 
c  * ( (IMOPT+IMAGIET)/IMAI) **2 . 

Mcap=ELECTEI/ 1000 . /specCap 

C  assume  lanode  and  cathode  are  tubes  with  thickness  20%  of  radius 
r eactor =Lanode*pi*rhomat* (Ranode**2 . +Rcath**2 . * ( 1 . 2**2 . -i . ) ) 

C  take  from  dimensions,  material  density,  etc 
injector=0. 

*****  EID  DPF  MASS  CALC  ***** 

BID  IF 

IF(reactype . EQ . 3)  THEM 
*****  BEGII  FRC  MASS  CALC  ***** 

theta-rhomat*Lcham*pi* ( (rcoil+drcoil) **2 . -rcoil**2 . ) 

C  [kg] 

plug=rhomat  *dzplug*pi* ( (rplug+drplug) **2 . -rplug**2 . ) 

C  [kg] 

reactor=theta+plug 

C  DPF  only  requires  puff  fill  gas  introduction  system,  not  injector 
C  assume  injector  mass  scales  linearly  with  mass  flow  to  be  injected 
in j  ector=mdf uel*specln j 
Mcap=0 . 

*****  EID  FRC  MASS  CALC  ***** 

EID  IF 

*****  EID  SPECIFIC  REACTOR  MASS  CALC  ***** 


*****  BEGII  GEIERIC  MASS  CALC  ***** 

mdout=mdfuel+mdprop 

ex=ElP(deltav/(Isp*g))-l 
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o  o 


t ms =ms  y s  *ex/ ( mdout * ( 1-Xtank*  ax ) ) 


nozzle=rhomat  *dznoz*pi*  (  (  rnoz+draoz )  *  *2 .  -rnoz*  *2 . ) 

C  [kg] 

as suae s  nozzle  has  rectangular  cross  section 
dznoz  *  drnoz 

C  shielding  Iron  parasitic  DDn 

C  assuaes  constant  flux  over  tms,  ending  with  dose=naxdose 
C  set  nflux  allowable  by  mu  dose 

nf lux=oaxdoae/ (tms*remnf lux) 

C  [cnT-2  s"-l] 

s igef f =rhoLiH/MVLiH* ( 1 6 . + 1 . ) *1 . e-24 
C  sigeff=[cm*-2]  effective  e  of  LiH  shield  material 

C  l.e-24  converts  bams  to  cm*2 

C  assume  nf lux  is  attenutated  as  nflux=nfluxO*EXP(-sigeff*drshld) 

C  solve  for  drshld 

nfluxO=RRddn*volplas*Xcone/(4. *pi*rpssngr**2. ) 

C  [s*-l] 

drshld=l .  /aigef  f  *L0G(nf  luO/nf  lnx) 

C  [cm] 

r2=rl>drshld 
C  [cm] 

volshld=4 . /3 . *pi* (r2**3 . -rl**3 . ) *Xcone 
shield=rhoLiH*volshld 
C  [kg] 

radiator=(l .  -w)*Pth*specRad 

C  assume  generator  mass  scales  linearly  with  power  required 
C  assume  2  specific  energies:  generators,  and  external  supply 
powersys=specGen*PelGen+specExt*PextReq+Mcap 


msys=injector+shield+radiator+powersys+Mcap+reactor+nozzle+ 

tpayload+shelter 


C  mass  of  fuel,  propellant  +  tanks  +  pumps,  etc. 
c  deltav  will  ultimately  come  from  mission  analysis  program 

fuel=tms*mdfuel 
prop=tns *mdprop 
tanks=Xtank* (f uel+prop) 
mf 3asys+tanks 
mO=mf+prop+fuel 

*****  EID  GEIERIC  MASS  CALC  ***** 

E1D 

C  subroutine  mass 
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A. 16  output. for 


C23466789  0 

SUBROUTINE  output 

INCLUDE  comblk.for 

************  OUTPUT  ****************** 

tanks=Xtank*(prop+fuel) 
tmd=tms/ (24 . *3600 . ) 

C  tmd [d] *3600 . *24 . =tms (s] 


C  display  output  values 
IF(loglin.EQ.l)THEN 
C  increment  IncVar  LOG 10 

MRITE(* ,92) 

k  ’+  heps' , count/howmany *100, ’%  ’ .IncVar, ’=' ,10**operate 

ELSE 

C  otherwise  increment  IncVar  linearly 

WRITE(*,92)  ’+  heps’ ,count/howmany*100, ’*/,  ’  .IncVar,  ’  =  ' , operate 

END  IF 


IF(ol)WRITE(55, 191) ’  T[keV]  =  ’,T 

IF(o2)THEI 

IFC.H0T.o9)  THEN 
WRITE(S6,*) 

WRITE(55,*) 'PLASMA  PARAMETERS’ 

IF(reactype . Eq . 2) THEN 

WRITE(56, 191) ’  Bnoz[T3=’ .Bnoz 
WRITEC55, 191) ’  pulse  voltageCV]*’ .volt 
VRITE(56 , 191) '  pulse  current [A] s', Imax 
ELSEIF (r eactype . Eq . 3 ) THEN 

HRITE(5B, 192) ’  Bnoz,  Bplug [T] =’, Bnoz Bplug 
VRITE(S6, 192) '  Bext,  Bint  CT]=’ .Bext , ’ , ’ ,Bint 
VRITE(66, 191) ’  beta  =’ .beta 
WRITE(B5 , 191) ’  taul  [s]  =’,tauN 
ENDIF 
END  IF 


IF  (Pinj.Eq.O)  THEN 
IF(o9)THEl 

VRITE(66, 131)T 
ELSE 

HRITE(B6, 191) ’  qp,rpPi/Pi=  1  *ignited*’ 
ENDIF 
ELSE 
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IF(oS)THES 

WRITE(S5, 132)T,TFP/Pinj 
ELSE 

WRITE(55 , 191) ’  Qp,rpPf/Pi=* ,TFP/Pinj 
EIDIF 
EIDIF 
BID  IF 


C  fuel  +  propellant  masses 
IF(o3)THEI 
IF(o9)THEI 

WRITE (55, 134) T.mdout, prop. fuel 
ELSE 

WRITE (55,*) 

WRITE (56 , * ) ' EXHAUST  MASSES  * 

WRITE(56, 191) ’  mdout[kg/s]=’ .mdout 
WRITEC55 , 191) ’  prop[kg]+sys  =’,prop 

WRITEC55, 191)  ’  fuel[kg]+sys  =’,fuel 

EIDIF 
EIDIF 


C 


component  masses 
IF(o4)THEI 
IF(o9)THEI 


WRITE ( 55 , 139) T , inj  ector , radiator , shield , reactor , nozzle , 


Rposersys .payload, shelter 


ELSE 

WRITE (55,*) 

WRITE(55,*) 'COMPOIEIT  MASSES’ 
WRITE(55, 191) ’  fuel  inj . [kg]  = 
WRITE(55 , 191) ’  tanks [kg] 
WRITE(55, 191) ’  radiator [kg]  = 
WRITE(55, 191) ’  shield [kg] 
WRITE(55, 191) ’  reactor[kg]  - 

WRITE(56 , 191) ’  nozzleCkg]  = 

WRITE(55, 191) ’  poser  sys.[kg]  = 

WRITE(56, 191) ’  payload [kg] 
WRITE(55 , 191) ’  shelter [kg]  = 

WRITE (56,*)  ’ - 


’  .injector 
' .tanks 
’ .radiator 
’ , shield 
’ .reactor 
’  .nozzle 
’ .posersys 
’ .payload 
' .shelter 


WRITE(65, 191) ’  final  oass  Dtg]=’,af 


EIDIF 

EIDIF 


C  output  posers  (heat) 

IF(o6)THEI 

IF(o9)THEI 

WRITE ( 55 , 134) T , Pradsall , Pcoil , Pnoz 
ELSE 

WRITE (55,*) 

WRITE (56,*) ’THERMAL  OUTPUT  POWERS  (EEAT) ’ 
WRITE(55 , 191) ’  Pradsall  [W] Pradsall 
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WRITE(55 , 191) ’  Preactor  [W]=', Pcoil 

WRITE (55, 191) '  Pnoz  [W]  =  ' .Pnoz 

WRITE (55,*)  * - * 

WRITE(55, 191)  ’  Pth  [W]=',Pth 

BID  IF 
E1DIF 

C  input  posers  (electrical) 

IF(o6)THEI 

IF(o9)THEI 

WRITE (55 , 136)T, Ppayld, Pin j/etalnj , Pcoil , Pplug ,Pnoz 
ELSE 

WRITE(65,*) 

WRITE(55,*) 'ELECTRICAL  ISPUT  POWERS' 

WRITE (55 , 191) '  Ppayld  [W]=', Ppayld 
WRITE(55 , 191) '  Pinj/oinj [W]=' .Pinj/etalnj 
WRITE(55, 191) '  Pcoil  [W]*', Pcoil 

WRITE(55 , 191)  ’  Pplug  [W]  =  \ Pplug 

WRITE(56, 191) '  Pnoz  [W]=',Pnoz 

WRITE (56,*)  ' - ' 

WRITE(S6, 191) '  PelReq  CW]=',PelReq 
EIDIF 
EIDIF 

C  other  posers  (electrical) 

IF(o7)THEI 

IF(.I0T.o9)THER 

WRITE(55,*) 

WRITE ( 55 . * ) ' GEIERATED  AID  EXTERSALLY  REQUIRED  ELEC.  POWER' 
WRITE(55,191) '  PelGen  CW]=' .PelGen 

WRITE(55 , 191) '  PextReq  [W]=» .PextReq 

WRITE(56, 191) '  ga*Gen=' .ganGen 
EIDIF 

IF  (PextReq. EQ.O)  THEM 
IF (o9) THEM 

WRITE ( 56, 131) T, ganGen 
ELSE 

WRITE(55, 191) '  Qp,spPelGen/PertReq=  1  *ignited*' 

EIDIF 

ELSE 

IF(o9)THEI 

WRITE (55 , 133)T .ganGen , PelGen/PextReq 
ELSE 

WRITE(55,191) ’  qp,spPelGen/PextReq=' , PelGen/PextReq 
EIDIF 
EIDIF 
EIDIF 

C  rocket  output  data 
IF(o8)THEI 
IF(o9)THEI 
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WRITE(S5, 134)  T.  thrust,  IFTHRST,  thrust /IFTHRST 
ELSE 

WRITE (56,*) 

WRITE(55,*) 'RAW  ROCKET  OUTPUT  DATA* 

WRITE(55, 192) »  TmixCK] ,Taix[eV]=* ,Taix, * , * ,Tmix/ 11804. 85 
WRITE(55, 191) *  n_aixEca'-3]=',n_aix 

WRITE(55, 191) *  Degree  of  propsllsnt  ionization  ** .ionized 
IF  (ionized.LT.0.3)  WRITB(S6,*)*  *  lot  ionizsd  enough  for  *, 
A 'magnetic  nozzle  ** 

IF  (Taix.LT.2000)  WRITE(55,*)'  *  use  noraal  nozzle  *' 

WRITE (55 , 192) '  Thrust Dl] .Thrust [lbf ]  =  * .thrust , * , ' .thrust/ 
A4.4482 

WRITE(55 , 191) ’  Specific  Impulse Cs]=' , lap 
WRITE(55, 191) ’  Mission  firing  time  Cd]  =  ',tmd 
WRITE (56, 192)'  MO.Mf  Ckg] ' ,a0, * , ' ,mf 
WRITE(56,*) 

WRITE(65,*) 'ROCKET  FIGURES  OF  MERIT' 

WRITE(55, 191) '  Initial  Thrust/weight  =* .thrust/ (aO*g) 
WRITE(56, 191) '  Jet  Power  w/o  fusion  [W]  =  \ 
AIFTHRST**2/(mdfuel*2) 

WRITE (56, 191) '  Jet  Power  w/  fusion  [W]=* ,g/2.*thruat*Iap 
WRITE (65 , 191) *  Fus.  Chrg.  Part.  Pow.  W]=',TFP 
WRITE(55, 191)  ’  Jet  Specific  Power  ‘0  [kg/kWjet]*' , 

Aa0/(g/2. *thrust*Isp/l . e3) 

WRITE (55 , 191) *  Overall  Efficiency  Pjet/Pf=* , 
k (g/2 . *thrust *Isp) /TFP 

WRITE(55 , 191) ’  Rocket  Eq.  Delta  Vee  Db/s]=* .deltav 
WRITE(56, 191) ’  Payload  aass/Systea  mass*' .pay load/aO 
WRITE (55 , 191) *  Mrpadprop/adf uel* ’ .adprop/adfuel 


EIDIF 

EIDIF 

IF(.I0T.o9)  THEM 
WRITE (55,*) 

WRITE(56,*) ******  lext  temperature  increment  ****** 
WRITE(65,*) 

EIDIF 

C  varitop  input  file 
c  IF(olO)THEI 

c  write(*,*) 

c  write(*,191)'  pO[kw, jet]** ,g/2.*thrust*Isp/l.e3 

c  write(*,191) *  mOCkg]=’,aO 

c  write(*,191) ’  isCs]**,Isp 

c  write(*,191) *  alfalCkg/ks]*’ ,af/(g/2.*thrust*Isp/l.e3) 

c  write(*,191) ’  alfa2[kg]*' ,aO-af 

c  write(*,191) '  tendCd]*’ ,tad 

c  write(65,191)’  pO*’ ,g/2.*thrust*Iap/l.e3 
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c  write(65,191)’  mO=',mO 

c  write(65, 191) ’  ia=',Isp 

c  write(6S,191)’  alial=’ ,mi/(g/2.*thrust*Isp/l.e3) 

c  write(65, 191) ’  alf a2= ' ,mO-ai 

c  urite(65, 191) ’  tend=',tmd 

c  EVOIF 

IF(ol2)VRITE(56 , 134) 

RT , del t  av , tad , payload/ mO 

c  Pprop=Pmix-Plp 

IF(ol3)WRITE(56 , 136) 

*T.IFTHRST**2/(mdluel*2) ,Plp-IFTHRST**2/(mdtuel*2) .Pmix-Plp, 
RPmix , g/2 . *thruet*Isp 

IF(ol4)URITE(56, 135) 

RT ,IFTHRST**2/ (mdfuel*2)/Pmix , (Plp-SFTHRST**2/(mdiuel*2) ) /Pmix , 
k (Pmix-Plp) /Pmix , g/2 . *thrust  *Isp 

ZF(olS)VRITE(55, 136) 

RT ,g/2 . *thrust*Isp , P j  etIF ,g/2 . *thrust*Isp/P j  etIF, 

JkPinj  ,g/2 .  *thrust* Isp/P in j 

E1D 

C  subroutine  output 
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A. 17  reactor.hep 


reactype:2=Perform  DPF  simulation, 3=Perform  FRC  simulation 

2. 

IncVar:  indsp. ,  incremented  var.:  T.deltav, . . 
dsltav 

howmany  =  points  to  graph  between  lower  and  upper  values  of  IncVar 
10 

lower  *  lower  value  of  IncVar 
5e2 

upper  =  upper  value  of  IncVar 
l.Be4 

loglin=  L0G10  or  linear  increment  of  IncVar  ( 1=L0G10 , other=linearly ) 

1. 

T  =  CkeV]  Plasma  thermal  energy  T=Ti+Ta 
160. 

sprot  =  fuel  particle  density  fraction  of  protons 
0. 

zdeut  = 

0.60 
ztrit  = 

0. 

x3He  = 

0.60 
z6Li  = 

0. 

xllB  = 

0. 

zalfa  = 

0. 

rt  =  unitless,  ratio  T(electron)/T(ion) 

1. 

Kc  =  unitless, fraction  of  cyclotron  radiation  absorbed  in  wall 

0.1 

cycabs  =  fraction  of  cyclotron  radiation  absorbed  in  plasma 
0.6 

y  *  unitless.frac.  charged  part,  power  absorbed  in  plasma 

0.8 

Xn  =  Ratio  of  fuel  density  to  propellant  density  in  reaction  chamber 
1. 

etaGen  =  generator  efficiency  (thermal->electrical) 

0.3 

rlCcm]  =  inner  radius  of  shield 

100. 

rpssngrCcm]  *  passenger  distance  from  source 

1000. 

Xcone  *  fraction  of  isotropically  radiated  neutrons  to  be  shielded  against 
0.126 

specRadCkg/UthJs  specific  mass  of  radiators  (620  kg/MVth) 

0.621e-3 

specGen[kg/Wal]=  specific  mass  of  generators 
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2. e-4 

specExt Ckg/Wel] =  specific  mass  of  external  electrical  eng.  supply 

3. e-4 

resmat [ohm  m]=resistivity  of  theta  pinch/mirror  materials  (Cu=1.72e-8) 

1.72e-8 

rhomsAt Ckg/cm'3] =density  of  theta  pinch/mirror  materials  (Al=2.67g/cm“3) 

8 . 9e-3 

tens max [Pa] = tensile  strength  of  theta  magnetic  materials(280e6sCu)  w/cables,etc. 
280.eS 

Xtenssratio  max  tens,  strenth  of  material  to  tensile  stress  (safety  factor) 

2. 

maxdose  [cSv]  over  mission  travel  time  (assuming  reactor  sole  source) 

S. 

ddnapin  =  suppressing/enhancing  factor  from  spinpolarization  ( l=no  effect) 

1. 

ddpspin  =  suppressing/enhancing  factor  from  spinpolarization 

1. 

d3Hespin  =  suppressing/ enhancing  factor  from  spinpolarization 

1. 
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A. 18  output. hep 


oO:WRITE  Initial  values  (T,t=yes,  F,l=no) 

F 

oi: WRITE  TiCkeV] 

F 

o2: WHITE  plasma  parameters 
F 

o3 : WRITE  exhaust 
F 

o4: WRITE  component  masses 
F 

06: WRITE  waste  output  posers  (heat) 

F 

06: WRITE  electrical  input  posers 
F 

o7: WRITE  generated  and  externally  required  elec,  poser 
F 

08: WRITE  Rocket  output  data 
F 

o9:WRITE  data  in  table  form  (without  text  explanations) 

T 

olO: WRITE  (was  VARITOP) 

F 

oil: 

F 

ol2:T,  deltaV,  tmd,  payload/aO 
T 

ol3 : T . mdfuel*Vrun**2/2 , Plp-adfuel*Vrun**2/2 , Pprop , Paix , P j  et 
F 

ol4 : T , P j  etIF/Paix , Pc/Pmix , Pprop/Paix , P j  et 
F 

ol5:T,Pjet,Pjet  s/o,  ratio 
F 

************************************************************ 

Key  for  table  form  output:  (incomplete) 
ol: 

o2:  T.  Qp 

o3:  T.mdout, prop, fuel 

o4 :  T , inj ector , radiator , shield , reactor , nozzle , posersys , payload , shelter 

06:  T,Pradsall,Pcoil,Pnoz 

06:  T.Ppayld.Pinj/etalnj .Pcoil.Pplug.Pnoz 

o7:  T,TFP/Piaj 

08:  T,  IFTHRST,  thrust,  thrust/KFTHRST ,  Isp, 
oll:T,  («FTHRST**2/(mdfuel*2)),  g/2.*thrust*Isp, 
ft (g/2 . *thrust*Isp)  /  (IFTHRST**2/(mdfuel*2)), 

*(g/2. *thrust*Isp)/TFP 
ol2:T,  deltaV,  tad,  payload/aO 
Plp*TDBLTAP+mdfuel*Vrun**2/2 

ol3 : T , mdf uel*Vrun**2/2 , Plp-adf uel*Vrun**2/2 , Pprop , Pnix , P j  et 
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o 14 : T , P j  atlF/Pinix , Pc/Pmix , Pprop/Pmix , P j  at 
ol5:T,Pjat,Pjat  s/o,  ratio 
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A. 19  dpf.hep 

RAVODE=AVODE  RADIUS  [cm] 

S.08 

RCATH=CATHODE  RADIUS  [cm] 

8.00 

LAIODE-AIODE  LEIGTH  [cm] 

38.2 

RH0I0V-IVITIAL  FILL  GAS  DEISITY  Ckg/ca‘3] 

2.2E-10 

V0LT=CAPACIT0R  BAIK  CHARGIIG  POTEXTIAL  IV  VOLTS 
27000. 

CAP-CAPACITAVCE  IV  FARADS 
3.6SE-4 

LIVIT=IVITIAL  CAPACITAVCE 
2.SE-8 

FSVPLV=SVOVPLOW  EFFICIEVCY  FACTOR 
0.7 

FPVCH=PIVCH  EFFICIEVCV  FACTOR 
0.25 

LPVCHsPIVCH  LEIGTH  [cm] 

2.54 

PVCHRADSPIVCH  RADIUS  [cm] 

1.5E-1 

REPRATE=PLASMA  FOCUS  FIRIVG  RATE  IV  S**-l 

100. 

PVCHTIM=DURATIOV  OF  STABLE  PIVCH  IV  SECOVDS 
l.OE-4 

ITERS=VUMBER  OF  ITERATIOVS 
500. 

DSCHRG=dur at ion  of  discharge  [a] 

1.0E-7 

IMAGVET=magnetic  currant  [A] 

3.18E5 
Bnoz* [T] 

8.2a3 

apacCap-CAPACITOR  BAIK  SPECIFIC  EVERGY  IV  KJ/KG 

2.0 

Xrnoz=ratio  nozzla  radiua/coil  radiua 

0.1 

apeclnj [kg/Wmdf ual] =  apacific  nut aa  of  fusion  fuel  in j actors (dpf=0, gas  puff) 

0. 

etalnj= Injector  efficiency  (electrical->inj acted) 

1. 
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A. 20  rocket. hep 

Nr  =  Ratio  of  propellant  mass  flow  to  fuel  mass  flow 
l.e3 

Zprop  =  unitless. atomic  number  of  propellant 

1. 

iprop  =  unitless .atomic  mass  of  propellant 

1. 

hO  »  [J/kg]  initial  enthalpy  of  propellant  (read  from  chart)  (SO  K) 

1 . 403e8 

hNat  3  [J/kg]  enthalpy  of  propellant  at  materials-limited  temperature  (2000K) 
1.716e8 

Xtemp  3  Ratio  of  mixing  temperature  (absolute)  to  throat  temperature  (abs) 
1.36 

Xvel  2  Ratio  of  exit  velocity  to  throat  velocity  of  magnetic  nozzle 

2. 

deltav  Cm/s]  =  mission  deltav 
3e3 

payload [kg] =  mass  of  desired  payload 
I.eS 

Ppayld[V]=  electrical  requirements  of  payload 
40.  e3 

shelter  =  [kg]  mass  of  Solar  Particle  Event  (SPE)  storm  shelter 
$.5e3 

dznoz= [cm] nozzle  thickness 

5. 

volmix= [cm* 3]  mixing  chamber  volume 
1 .  e7 

Xtank=tank  fraction,  determine  mass  of  tanks,  feed  lines,  etc. 

0.1 

Tmix [K] =initial  guess  for  exhaust  mixing  temperature 
7.  e3 
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A.21  tm.for 


program  tm 


implicit  real  (a-z) 

CHARACTER  pause* 1 

DATA  pi , Wd3He , grav/3 . 141S926 ,18.3,9.8/ 

131  FORMAT  (1PE11.3) 

132  FORMAT  (1P2(E11.3)) 

133  FORMAT  (1P3(E11.3)) 

134  FORMAT  (1P4(E11.3)) 

135  FORMAT  (1P6(E11.3)) 

137  FORMAT  (1P7(E11.3)> 

138  FORMAT  (1P8(E11.3)) 

141  FORMAT  (A.1PE11.3) 

142  FORMAT  (A, 1P2(E11 . 3) ) 

143  FORMAT  (A, 1P3(E11 . 3) ) 

144  FORMAT  (A,1P4(E11.3)) 

145  FORMAT  (A,1P5(E11.3)) 

147  FORMAT  (A,1P7(E11.3)) 

148  FORMAT  (A, 1P8(E11 .3)) 

Pi(kT)  ■  inven*inven/Vp/4 . *Wd3He*SVd3He (kT) *1 . 602ell  !MJ/us 

Pbr(kT)  »  2.5*inven*inven/Vp*SQRT(kT)*Se-7  !MJ/us 

g(kT)  =  (Um*le-3*Im  +  P*(kT)  -  Pbr(kT))/(inven*1.602e-4)  !keV/us 

OPEI  (UHIT=46,  FILE=’tmi.dat' ,  STATUS3 ’ OLD ' ) 

READ (45,*) 

READ(46,*)W 

READ(46,*) 

READ (46,*) Ido tO 
READ  (45,  *) 

READ(46,*)Ldot 
READ(46, *) 

READ (46 , * ) inven 
READ  (45,*) 

READ(45,*)rc,ra 
READ(46, *) 

READ(46,*)nu 
READ (45,*) 

READ ( 45 , * ) rt ank 
READ(45,*) 

READ(46,*)mpay 
READ (45,*) 

READ(45,*)deltav 
CLOSE  (U1IT=45) 

urit*(*,141) ’  V  (MJ)',W 

»rite(*,141) *  IdotO  (MA/us) ' ,IdotO 
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write(*,141) ' 
write(*, 141)  ’ 
write(*, 142) ’ 
write(*,141)  ’ 
write(*, 141) ’ 
write(*,141) ’ 
write(*,141) ' 


Ldot  (aH/ a ) ’ , Ldot 

inventory  I  (*lel8) ’ ,inven 


rc.ra 
reprate  nu 
ztank 
apay 
deltav 

WRITE(*,*) 'Press  <Enter>  to  continue' 
READ(*, ’ (Al) ')pause 


(ca) ' , rc.ra 
(Hz) ’ ,nu 
' , ztank 
(kg)', apay 
(ka/s) ' , deltav 


speccap=2e-3 

acap=W/speccap 

write(*,141)  ’  acap  (kg)'.acap 

AU= ( W*IdotO*Ldot*Ldot )**(l./3.)*lel 

AI=(W*IdotO/Ldot)**(l./3.)*iel 

At= ( V/ ( IdotO* ldot 0*Ldot ) ) ** ( 1 . /3 . ) * 1 e 1 

write(*,143) '  AU(kV) ,AI(HA) ,At(us) ' .AU.AI.At 

Ua=0 . 64* AU 

la=0 . 64* AI 

write(*,144) '  UO.Ua(kV) .Ia(MA) .trise(us) ' ,2. 12eAU.0a.Ia, 1.43*At 
tau_p=0 . 268* At 

write(*, 141) ’  tau_p  (us)',tau_p 
lp=rc-ra  ! cm 

kTO=Im*Im*lp/ (2*inven*l . 602e-l) 

write(*, 141) '  kTO  (keV)  ' ,kT0 

write(*, 141) ’  SV(kT0)(cc/s) ' ,SVd3He(kT0) 

kr=70. 

rO=ra/kr 

Vp=pi*rO*rO*lp  !ca"3 
dens=inven*lel8/Vp  !ca"-3 

write(* , 143) ’  rO(ca) ,Vp(cc) ,dens(/cc) ' .rO.Vp.dens 

write(*, 143) '  Umla.Pl.Pbr  (MJ/us) ' .Ua*le-3*Ia,Pf (kTO) .Pbr(kTO) 

Pis  v= inven*inven/ Vp/4 . *tfd3He* 1 . 602e 1 1 

PbrkT=2 . S* inven*inven/Vp*6e-7 

write(*, 141) '  Pl(kTO)/<sv>  =’,Plev 

write(*, 141) ’  Pl(kTO)/<sv>/I  =' ,Plsv/(inven*1.602e-4) 

write(*, 141) '  Pbr(kT0)/kT*l/2  =',PbrkT 

write(*,141) ’  Pbr(kT0)/kT'l/2/I  =' ,PbrkT/(inven*1.602e-4) 

h=tau_p  !us 

gl=g(kTO) 

g2=g(kT0+h/2 . *gl) 

g3=g(kT0+h/2. *g2) 

g4=g(kT0+h*g3) 

write (* , 143) ' kTl , kT2 ,kT3 ' ,kTO  fa/2 . *gl , kTO+h/2 . *g2 ,kT0+h*g3 
grk  *  (gl+2 . *g2+2 . *g3+g4)/8 . 
kTl  =  kTO  +  grk*b 

write(*,145) '  gl,g2,g3,g4,grk’ ,gl,g2,g3,g4,grk 

kTave  »  kTO  +  grk*h,2. 

write (*,142) '  kTl, kTave  (keV)' , kTl .kTave 

El  =  tau_p*Pl (kTave)  !HJ 

Ebr  =  tau_p*Pbr (kTave)  !MJ 
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write(*, 144) ’  Ef.Ebr  (MJ) ,  Qf  ,qbr\ Ef.Ebr  .Ef/V.Ebr/V 
Hf_H  =  (2  *  Ef)/(Wd3He*inven*1.602e-l) 
write(*, 141) *  Hf/H  *,Sf_S 

IF  (Hf_H.GT. 1 .0)  write(*,141) '  Qmax=’ , inven*Wd3he/W*le-l 

write(*,*) 'au?' 

read(*,*)nu 

Pth  =  («  +  Ef  +  Ebr)*nu  !MW 
dal tab  =  1.05e2  !MJ/kg 

c  deltah  =  31.3  !MJ/kg  using  old  1.716-1.403  e8J/kg 

mdot  =  Pth/deltah  !kg/s 

srite(*. 142) ’  Pth(MW) ,adot(kg/s) * , Pth, mdot 
Fprop  =  Pth  -  W*nu/0 . 5  !MW 
Vax  =  SQRT (2*Pprop/ mdot )  !ka/s 
write(*,142) *  PpropOW) ,Vex(km/s) ' ,Pprop,Vex 
kTprop=Prop/mdot*l . 67a-2/l . 602 

write(*. 141) ’  kTprop.kTthr  (eV) ’ ,kTprop,kTprop/1.3S 
thrust=adot*Vex 

write(*, 142) ’  F(kH) ,Isp(*10*3  s) ’ .thrust ,Vex/grav 
write(*,*) 'Deltav.aax  (ka/s) 1 ,vex*LOG(l+l./xtank) 
c  write(*,*) 'Delta  v?’ 

c  read(*,*)deltav 

ms=acap+mpay 
oe=EXP(deltav/vex)  -  1 
tmis8=ms/mdot*ee/(l-xtank*«e)  !s 
tad=taiss/8 . 640e4 
writa(*, 141) ’  tmias  (d),,tBd 
mO=BB+mdot*taiss*(l+xtank) 
alpha=mO/ ( Pprop* 1 e3 ) 

write(*, 143) '  aO(kg) ,npay/aO, alpha (kg/kV) ’ .aO.apay/aO, alpha 

f_w=thrust*le3/(a0*grav} 

write(*,141) ’  F/«0 ’  ,f_w 

OPES  (UIIT=46,  FILE=’tao.dat' ,  STATUS =' OLD ' ) 
write (46 , 137)apay , deltav , nu , tad, f _w , alpha , apay/aO 
end 


C  for  3He(d,p)‘  rxn;  ca*3/s 
REAL  FUSCTIOI  SVd3He(T) 

REAL  a.b.c,d,T.x.L0G10 
x=L0G10(T) 
a=  0 . 353S689 
b=  -3.310354 
e*  10.10471 
d=-25. 67344 

SVd3He=10**(a*x**3+b*x**2+c*x+d) 

ESD 

C  function  SVd3He 
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A. 22  tmi.dat 


W  (MJ)  =  ?M 
100 

IdotO  (NA/us)  =  ?" 
10 

Ldot  (mH/a)  =  ?" 

1 

invan  (*lal8)  =  ?" 
100 

rc,ra  (cm)  =  ?" 

8  3 

rapr&te  nu  (s"-l)" 
1 

xtank=?" 

.1 

mpayload(kg) 

leS 

daltav  (km/ a) 

28 


B  c  source  code 


B.l  cir.c 


♦include<math . h> 

#include<atdio . h> 

♦define  muO  1.2566E-6 
♦define  pi  3.141S926S4 
♦define  J.MeV  1.60219E-13 
♦define  J.keV  1.60219E-16 
♦define  Vddn  3.27 
♦define  Uddp  4.03 
♦define  Wd3He  18.3 

FILE  *fp,*fq; 

float  time[1000] ,cv[1000] ,i[1000] ,find_max(); 
int  max.records; 

void  graph () ,read_data() , enter. variables () , across ; 
char  xtype [20] , ytype [20] , ans ; 
float  x.data [1000] ,y_data [1000] ; 

float  A[ll];  /a  polynomial  solution  coefficients  */ 

/a  C  arrays  start  at  0  in  references,  bat  not  in  definitions  */ 
void  enter. variables () ,description() , rundown () ; 
float  UstarO  ,Istar()  ,SVd3He() ; 
float  tO.tf .tdivs.U.IdotO.Ldot; 

float  U0 , t , h , Uopt , Imopt , t opt , tmax , Imax , uuopt , Umax , L0 , Lmax , 

Cap , OmegaO , nopt , lambda , Rtot , tan , Omega ,  Qmax , Bcoef f , TFE , TEBrem , 
deut , hel3 , RRddn , RRddp , RRd3He , dPFddn .dPFddp , dPFd3He , dEf tot , dEbrem , 
volplas , kT.keV ; 
int  test; 
iat  j; 

float  Wrundown , W.pinch , t.pinch ; 

float  kT, inventory ,Zeff ,r_0,r_anode ,r_ cathode ,r_p inch, 

compr ession.rat io , vol.pinch , dens it y , P.brem , P_f us ion , 

E_br em , E.f us ion , Q ; 

float  AU , AI , At , t star , UOstar , Imstar , Umstar , tmstar ; 
float  current , voltage ; 


mainO 

{ 

int  first_time=l; 
do 

{ 

clrscrQ ; 

printfC'DPF  Circuit  Program\n\n") ; 
if  (first.time)  descriptionO ; 
first.time=0; 
rundown 0 ; 

printf ("\n\nD0  YOU  WISH  TO  ALTER  DATA  AID  C0MHEICE  AI0THER  RUI  <Y/I>  "); 

> 
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while(ana=getch()==’y '  II  ana=='Y’) ; 
clracrO ; 

>  /*  main  */ 


void  rundown () 

{ 

float  t«mp; 
float  vrd; 

enter_variables() ; 

vrd=Ldot*2*pi/ (muO*log(r_cathode/r_anode) ) ; 

AU=pow( (W*IdotO*Ldot*Ldot) , ( 1 . /3 . ) ) ; 

AI*pow( (W*IdotO/Ldot) , (1 ./3. )) ; 

At =pow ( (W/ ( IdotO«IdotO*Ldot ) ) , ( 1 . /3 . ) ) ; 

/* 

Uopt  a  2.12  *  AU; 

Imopt  =  0.64  *  AI; 
topt  =  1.6  *  At; 

*/ 

U0star=2. 12; 

U0=AU*U0atar; 

/*  GET  COEFFICIEITS  TO  SERIES  S0LUTI0*  */ 

/*  U*=SUH  A[j]t*'j  */ 

A [0] aUOstar ; 

A[1]=0. ; 

for  (j=0; j<=8; j++) 

{ 

ACj+2]  =  -ACj]*U0»tar/(2*(j+l)*(j+2)) 

-A[j+l]*(j+l)/<( j+2)*U0etar) ; 

> 

if ( ( f p=f open ( " c ir_out . dat " , "wt " ) ) **HJLL) 

{ 

printf ("Uboa!  Cannot,  open  output  file:  cir_out.dat\n") 
exit(l) ; 

> 

/*  find  maximum  current,  time  at  max  current  */ 

/*  RUIDOV*  PHASE  */ 

lmax»0; 

tmaxaO; 

/*  define  timestep  */ 
ha(tf-t0)/(tdive-l) ; 

/*  use  while  loop  here  */ 

VrundownaO ; 
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tstar=tO/At ; 


while  (Ustar(tatar)>0)  /*  quit  if  voltage  dropa  below  0  */ 

< 

voltage=AUsUstar(tstar) ; 
current=AI*Iatar(tatar) ; 
t*At*tstar; 
if  (current>Imaz) 

{ 

/*  get  the  initial  conditions  for  pinch  phase  */ 
tmax=t ; 

I*ax*current ; 

Umar=voltage ; 

Wrundown+=volt agescurrentsh ; 

if  (t>tf)  printf ("feed  longer  rundown  tine\nM); 

> 

f printf (fp, '7.1 1.3e  %11.3e  %11.3e\n",t, voltage, current); 
tstar+=h/At ; 

> 

f  dose(fp) ; 

if  (Wrundown  <=  V) 

{ 

W_pinch=W-Wnindown ; 

> 

else  { 

printf ("lo  current  maximum:  taking  V_pinch=0.1V\n") ; 
W_pinch=0.1*W; 

Imax=AI*Istar(0.9*tmax) ; 

Uaax=AU*Uatar (0 . 9*tmax) ; 

> 

t .pinch  =  tf_pinch/(Imax*Umax) ; 

kT=W_pinch/ ( inventory *J_keV) ;  /*keV*/ 

r_0=r_anode ; 

r_pinch=r_0/co*presaion_ratio ;  /stake  r0=ra*/ 

vol_pinch=pi*r_pinch*r_pinch*(r_cathode-r_anode) ; 

/stake  lpinch=rc-ra,  experimental  scalings/ 
dens ity=inventory/vol_pinch ; 
P_bre«=density*5e-37sZeffsdensity*sqrt(kT) ; 

P_fusion=denaity/4. sSVd3He(kT)*density*le-6sWd3Hesj_MeV; 
B_brea=»P_brea*vol_pinch*t_pinch ; 

E_fusion=P_fusionsvol_pinchst_pinch; 

Q=E_fusion/W; 

printf ("\n") ; 

printf ("vrd  (m/s)%11.3e  \nM,vrd); 

printf ("l.anode  (m)%11.3e  \n",vrd*tmax) ; 

printf  ("00s, Isa, ts«  y.ll.3e%11.3e%11.3e  \n",UO/AU,Imax/AI,tmax/At) ; 

printf ("taax.tf  (s)%il.3e%11.3e  \n",taax,tf ) ; 

printf ("Imax  (A)%11.3e  \nM,Iaax); 
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pr inti ( "UO , Umax 
print! ( "W , W rundown 
print! (MWrd/W 
print! ( "t_pinch 
print! ("\n") ; 
print! ("kT 
print!("P«ak  Brea, 
print! ( "r_0 , r .pinch 
print! ("density 
print! ("vol.pinch 
print! ("<ev> 

print! ("P_brea,P_lusion  (y)Xll.3e%11.3e\n",P_brea,P_lusion) ; 
print! ("E_brea,E_ins ion  (W ) y.l 1 . 3e%l 1 . 3e\n" , E_br ea , E_! us ion) ; 
print! ( "q , E_br ea/y  %1 i . 3e%l 1 . 3e\n" , q , E_br ea/y ) ; 


(V)'/.11.3eXll.3e  \n"  .UO.Uaax) ; 

( J)’/,l  1 . 3e'/,l  1 . 3e  W.y.yrundown); 

%11.3e  \n",yrundown/y) ; 

(s)%11.3e  \n",t_pinch) ; 

(keV)%11.3e\n".kT); 

(i)Xll.3e\n",kT/6.20); 

(a)%l 1 . 3e , %1 1 . 3e  \n" ,  r_0 ,  r .pinch) ; 
(a“-3)Xll.3e  \n", density) ; 

(m*3)y,11.3e  \n",vol_pinch) ; 

(  ca‘3/s )%1 1 . 3e\n" . SVd3He (kT) ) ; 


il ( ( ! q=iopen ( " c ir_ io . dat " , " wt " ) )==IULL ) 

{ 

print! ("Vhoa!  Cannot  open  output  !ile:  cir_io.dat\n") ; 
exit(l) ; 

> 

! print!  (!q,  '7.10 .  3e’/,10 . 3e'/.10 . 3e%10 . 3eXlO .  3e%10 . 3ey.l0 . 3eXlO .  3e" , 

V , IdotO ,Ldot , r.cathode , r.anode , inventory , Zell , coapression.ratio ) ; 
f  print!  (!q ,  ‘7.10 .3e%10.3ey.10. 3e%10 . 3e%10 . 3e%10 . 3e%10 . 3e%10 . 3e\n" , 
Iaax . t.pinch , r_pinch , kT , density , E_!usion . E.brea , q ) ; 

!close(!q) ; 


print! ("\nD0  YOU  VISE  TO  GRAPH  DATA  <Y/I>") ; 
i!  (ans=getch()=='y *  I |  ans=='Y’) 

read_data( ) ; 
test=menu() ; 
while (teat !*0) 

{ 

graphO; 
test=aenu() ; 

> 

> 

>  /+  rundown  */ 
void  descriptionO 

print! ("prograa  calculates  capacitor  voltage,  current  vs.  tine  \n") 

print! ("lor  siaplilied  rundown  phase  lor  DPF  circuit.  \n") 

print! ("capacitor  charge  *  series  solution  to  RLC  circuit,  L*L0+Ldot*t  \n") 
print! ("q«S0M  Ai*t**i  \n") 

print! ("10th  order  polynoaial  \n") 

print! ("prograa  author:  R.  Iachtrieb,  lov91  \n") 

print! ("last  aodilied:  Har92  \n\n”) 

} 
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float  Ustar(t) 
float  t; 

{ 

float  voltage; 

voltage  =  A CO]  +  A[l]*t  +  A[2]*t*t  +  A[3]*t*t*t  + 

AC4] *pos(t ,4)  +  A[S]*pov(t,6)  +  A[6]*poe(t,6)  +  A[7]*pou(t,7)  + 
AC8]»pOH(t,8)  +  A[9j*pov(t.9)  +  AClO>poe(t ,  10) ; 
return  (voltaga); 

> 

float  Iatar(t) 
float  t; 

float  current; 

current  =  -2./(U0star*U0star)* 

(A[l]  +  2*A[2]*t  +  3*A[3]*t*t  +  4*AC4] *t*t*t  +  5*AC5]*pow(t,4)  + 
6*A[6]*pow(t,S)  +  7*A[7]*pow(t,6)  +  8*A[8]*poe(t,7)  +  9*A[9] *poe(t,8)  + 
10*A[10] *pou(t ,9)) ; 
return  (current) ; 

> 


